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ABSTRACT 

We describe a photon-conserving radiative transfer algorithm, using a spatiahy- 
adaptive ray tracing scheme, and its parahel implementation into the adaptive mesh 
refinement (AMR) cosmological hydrodynamics code, enzo. By coupling the solver 
with the energy equation and non-equilibrium chemistry network, our radiation hy- 
drodynamics framework can be utilised to study a broad range of astrophysical prob- 
lems, such as stellar and black hole (BH) feedback. Inaccuracies can arise from large 
timesteps and poor sampling, therefore we devised an adaptive time-stepping scheme 
and a fast approximation of the optically- thin radiation field with multiple sources. 
We test the metho d with seve r al rad i ative transfer and radiation hydrodynamics tests 
that are given in llliev et aD (|2QQ6L |2Q09[ ). We further test our method with more 
dynamical situations, for example, the propagation of an ionisation front through a 
Ray leigh- Taylor instability, time-varying luminosities, and collimated radiation. The 
test suite also includes an expanding H ll region in a magnetised medium, utilising 
the newly implemented magnetohydrodynamics module in ENZO. This method linearly 
scales with the number of point sources and number of grid cells. Our implementation 
is scalable to 512 processors on distributed memory machines and can include radi- 
ation pressure and secondary ionisations from X-ray radiation. It is included in the 
newest public release of enzo. 
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> 1 INTRODUCTION 



Radiation from stars and black holes strongly affects their 
surroundings and plays a crucial role in topics such as 
stellar atmospheres, the interstellar medium (ISM), star 
formation, galaxy formation, supernovae (SNe) and cos- 
mology. It is a well-studied problem (e.g. [M athews 1965; 
Rvbicki k LightmanI IToTOI : iMihalas Mihalasa.1984 : .Yorkci 



methods: numerical — hydrodynamics — radiative trans- 



static uniform neutral medium. When recombinations bal- 
ance photo-ionisations, the radius of a so-called H ii region, 



Rs 



3N^ 



1/3 



(1) 



19861 ): however, its treatment in multi-dimensional calcula- 



tions is difficult because of the dependence on seven variables 
- three spatial, two angular, frequency, and time. The non- 
local nature of the thermal and hydrodynamical response 
to radiation sources further adds to the difficulty. Depend- 
ing on the problem of interest some simplifying assumptions 
may be made. 

An important case was considered bv lStromgrenI (|l939l l 
for an ultraviolet (UV) radiation source photo-ionising a 
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where is the ionising photon luminosity, is the re- 
combination rate, and uh is the ambient hydrogen num- 
ber density. Furthermore he found that the delineation be- 
tween the neutral and ionised medium to be approximately 
the mean free path of the ion ising; radiation . His semina l 
wor k was expanded upon by ISpitzed (|l948l . Il949l . Il954h 
and 'Spitzer Sav edofl (|l950 l). who showed that the ion- 
ising radiation heated the medium to T ^ 10^ K. If the 
density is equal on both sides of the ionisation front, then 
this over-press urised reg;ion would expand and drive a shock 
outwards fe.g. lOortlfTiH^ : ISchatzman Kahnlll955h . These 
early works provided the basis for the modern topic of ra- 
diation hydrodynamics of the ISM. A decade later, the first 
radiation hydrodynamical numerical models of H ii regions 
in spherical symmetry and plane-parallel ionisation fronts 
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were developed (e.g. lMathewslll965l : lLaskerlll966l : [Hiellminj 
1 19661 ). They described the expansion of the ionisation front 
and the evolution of its associated shock wave that carries 
most of the gas away from the source. At the same time, 
theoretical models of ionisation fronts matured and were 
classified by Kahn (1954) and Axford (11961 1) as either R- 
type (rare) or D-type (dense). In R-type fronts, the ionised 
gas density is higher than the neutral gas density, and in 
D-type fronts, the opposite is true. R-type fronts travel su- 
personically with respect to the neutral gas, whereas D-type 
fronts are subsonic. Furthermore "weak" and "strong" R- 
type fronts move supersonically and subsonically with re- 
spect to the ionised gas, respectively. The same terminology 
conversely applies to D-type fronts. "Critical" fronts are de- 
fined as moving exactly at the sound speed. These works 
established the evolutionary track of an expanding H ii illu- 
minated by a massive star in a uniform medium: 

(i) Weak R-type: When the star (gradually) starts to 
shine, the ionisation front will move supersonically through 
the ambient medium. The gas is heated and ionised, but 
otherwise left undisturbed. This stage continues until r ^ 
0.02Rs. 

(ii) Critical R-type: As the ionisation front moves out- 
wards, it begins to slow because of the geometric dilution 
of the radiation. It becomes a critical R-type front, which is 
equivalent to an isothermal shock in the neutral gas. 

(iii) Strong and weak D-type: The front continues to slow, 
becoming a strong D-type front, and then a critical D-type 
front. From this point forward, it is moving subsonically 
with respect to the ionised gas, i.e. a weak D-type front. 
Thus sound waves can travel across the ionisation front and 
form a shock. The ionisation front detaches from the shock, 
putting the shock ahead of the ionisation front. 

(iv) Expansion phase: After the shock forms, the H ii re- 
gion starts to expand, lowering the interior density and thus 
the recombination rates. This increases the number of pho- 
tons available for ionising the gas. The sphere expands until 
it reaches pressure equilibrium with the ambient medium at 

In the 1970's and 1980's, algorithmic and computa- 
tional advances allowed numerical models to be expanded 
to two dimension s, mainly using axi-symm e tric to simplify 
the problem (e.g; . Bodenheim er et all 1 19791 : ISand ford et al.l 
Il982l : lYorke et al. 1983). One topic that was studied exten- 
sively were champagne flows. Here the source is embedded 
in an overdense region, and the H ii region escapes from 
this region in one direction. The interface between the am- 
bient and dense medium was usually set up to be a constant 
pressure boundary. When the ionisation front passes this 
boundary, the dense, ionised gas is orders of magnitude out 
of pressure equilibrium as the temperatures on both sides of 
initial boundary are within a factor of a few. In response, 
the gas is accelerated outwards in this direction, creating a 
fan-shaped outflow. 

Only in the past 15 years, computational resources have 
become large enough, along with further algorithmic ad- 
vances, to cope with the requirements of three-dimensional 
calculations. There are two popular methods to solve the 
radiative transfer equation in three-dimensions: 

• Moment methods: The angular moments of the radi- 



ation held can describe its angular structure, which are 
relate d to the energy ener g y, flux, and radiati on pres- 
sure (jAuer MihalasI I197QI : iNorman et all Il998l ). These 
have been implemented in c onjunction with short char- 
acteristics ("stone et al." 1992, 2D), with long characteris- 
tics (Finlato^et al. 2009), with a variable Eddington ten- 
sor in the optically thin limit (OTVET; ICnedin Abell 
1200 ll : IPetkova SDringell[20Q9|) . and with an Ml closu re re- 



lation jGon5lez^ral]^^ Mo- 
ment methods have the advantage of being fast and indepen- 
dent of the number of radiation sources. However, they are 
diffusive and result in incorrect shadows in some situations. 
• Ray tracing: Radiation can be propagated along 



rays that extend throuj 


yh the computational grid (e.g. 


iRazoumov k ScottI Il999 


: Abe] 


et al.1 I1999I: Ciardi et al.1 


'2001': ISokasian et al.1 l200ll: 


Whalen & NormanI 20061: 



2006a: Trac k CenI 20071: iKrumholz et al.1 


I2OO7I: 


Paardekooper et al. [20101) or particle set (e.g 


. ISusa 


20061: 1 Johnson et al.1 12007: IPawl 


Lk k Schavd I2OO8 


, [2oTi 


Altav et al.l 20081: Hasegawa et al. 


I2OO9!). In general, these 



methods are very accurate but computationally expensive 
because the radiation fleld must be well sampled by the 
rays with respect to the spatial resolution of the grid or 
particles. 



Until the mid-2000's the vast majority of the three- 
dimensional calculations were performed with static den- 
sity distributions. One example is calculating cosmological 
reionisation by pos t-processing of density fields from N- 
body simulatio ns (Ciar di et al.ll 200l': 'Sokasian et al.l l200ll : 
iMcQuinn et~ah 2007: I liev et al.l l2006. 2007). Any hydrody- 
namical response to the radiation field was thus ignored. 
Several radiative transfer codes were comp a red i n four 
purely radiative transfer tests in llliev et aU (|2006l . here- 
after RT06). Only recently has the radiative transfer equa- 
tion been coupled to the hy drodynamics in three-dimensions 
e.g. Krumholz et al.|[2007l ) . In the second comparison paper 
Iliev et al .II2OO9'. hereafter RT09), results from these radia- 
tion hydrodynamics codes were compared. Even more rare 



are ones that couple it with magneto- hydrodynamics (e.g. 
IKrumholz et al.ll2007l ). The tests in RT06 and RT09 were 
kept relatively simple to ease the comparison. 

In this paper, we present our implementat ion, 
ENZO+MORAY, of adaptive ray tracing (Abel k Wandeltl 
l2002h in the cosmological hydro dynamics adaptive mesh 
refinement (AMR) code, enzo (jBrvan k NormanI Il997l : 
O'Shea et al. 2004). The radiation field is coupled to the 
hydrodynamics solver at small time-scales, enabling it to 
study radiation hydrodynamical problems. We have used 
this code to investigate the growth of an H II region fr om a 
IOOM0 Population III (Pop III) star dAbel et al.ll2007l). the 
early stages of reionisation from Pop III stars ( Wise k Abell 
2008a), radiativ e feedback on the formation of high red shift 
dwarf galaxies (jWise k Abell l2008bl : IWise et al.1 120101 ). ul- 
traviolet radiatio n escape fractions from dwarf galaxies be- 
fore reionisation (Wise k Cen 2009), negative radia tive feed- 
back from accreting Pop III seed black holes (Alv arez et al.l 
I2OO9I ). and radiative feedba ck in accreting supermassive 
black holes (|Kim et al.l I2OIIL in prep.). We have included 
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ENZO+MORAY in the latest public release of enzcQ, and 
it is also coupled with t he newly added MHD solver in 
ENZO (|Wang Abelll200^ y 

We have structured this paper as follows. In Section 
2, we describe the mathematical connections between adap- 
tive ray tracing and the radiative transfer equation. Further- 
more, we detail how physics other than photo-ionisation and 
photo-heating are included. We then derive a geometric cor- 
rection factor to any ray tracing method to improve accu- 
racy. We end the section by describing a new computational 
technique to approximate an optically-thin radiation field 
with ray tracing and multiple sources. In Section 3, we cover 
the details of our radiation hydrodynamics implementation 
in ENZO, specifically (1) the ray tracing algorithms, (2) cou- 
pling with the hydrodynamics solver, (3) several methods to 
calculate the radiative transfer timestep, and (4) our par- 
allelisation strategy. We present our results from the RT06 
radiative transfer tests in Section 4. Afterwards in Section 5, 
we show the results from the RT09 radiation hydrodynamics 
tests. In Section 6, we expand on these tests to include more 
dynamical and complex setups to demonstrate the flexibility 
and high fidelity of enzo+moray. Section 7 gives the results 
from spatial, angular, frequency, and temporal resolution 
tests. In Section 8, we illustrate the improvements from the 
geometric correction factor and our optically-thin approx- 
imation. We also show the effects of X-ray radiation and 
radiation pressure in this section. Finally in Section 9, we 
demonstrate the parallel scalability of enzo+moray. Last 
Section 10 summarises our method and results. 



2 TREATMENT OF RADIATIVE TRANSFER 

Radiation transport is a well-studied topic, and we be- 
gin by describing our approach in solving the radia- 
tive transfer equation, which in comoving coordinates 
(|Gnedin Ostriker|[l997l ) is 



1 dU n-Vh 



dt 



+ ■ 



(2) 



Here Ii, = /(z/, x, Q,t) is the radiation specific intensity in 
units of energy per time t per solid angle per unit area per 
frequency u. H = d/a is the Hubble constant, where a is the 
scale factor, d = a/aem is the ratio of scale factors at the 
current time and time of emission. The second term repre- 
sents the propagation of radiation, where the factor 1/d ac- 
counts for cosmic expansion. The third term describes both 
the cosmological redshift and dilution of radiation. On the 
right hand side, the first term considers the absorption coef- 
ficient = Hiu{y^^ t), and the second term ji, = juij^^ z^, t) 
is the emission coefficient that includes any point sources of 
radiation or diffuse radiation. We neglect any {v/c) terms 
in equation (|2]) that become important in the dynamic dif- 
fusion limit (lHi(v/c) ^ 1), where / is the characteristic size 
of the system. This occurs in relativist ic flows or very op- 
tically thick systems, such as stellar interior s or radiation- 
dominated shocks fsee iKrumholz et"aL I I2OO7I . for a rigorous 
derivation that includes (v/c) terms to second-order). 

Solving this equation is difficult because of its high 
dimensionality; however, we can make some appropriate 

^ http://enzo.googlecode.com 
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approximations to reduce its complexity in order to in- 
clude radiation transport in numerical calculations. Typi- 
cally timesteps in dynamic calculations are small enough 
so that Aa/a <^ 1, therefore a = 1 in any given timestep, 
reducing the second term to hdlu/dyi. To determine the im- 
portance of the third term, we evaluate the ratio of the third 
term to the second term. This is ifL/c, where L is the sim- 
ulation box length. If this ratio is ^ 1, we can ignore the 
third term. For example at z = 5, this ratio is 0.1 when 
L = c/H{z = 5) = 53 proper Mpc. In large boxes where the 
light crossing time is comparable to the Hubble time, then it 
becomes important to consider cosmological redshifting and 
dilution of the radiation. Thus equation (|2]) reduces to the 
non-cosmological form in this local approximation. 



1 diu . diu 
c at ax 



(3) 



We choose to represent the source term ji, as point sources 
of radiation (e.g. stars, quasars) that emit radial rays that 
are propagated along the direction h. Next we describe this 
discretisation and its contribution to the radiation field. 



2.1 Adaptive Ray Tracing 

Ray tracing is an accurate method to propagate radiation 
from point sources on a computational grid, given that there 
are a sufficient number of rays passing through each cell. 
Along a ray, the radiative transfer equation reads 

IdP dP 
c at or 



where P is the photon number flux along the ray. To sample 
the radiation field at large radii, ray tracing requires at least 
Nray = 4:7T / (Ax)"^ rays to sample each cell with one ray, 
where R is the radius from the source to the cell and Ax is 
the cell width. If one were to trace Nray rays out to R, the 
cells at a smaller radius r would be sampled by, on average, 
(r/R)^ rays, which is computationally wasteful because only 
a few rays per cell, as we will show later, provides an accurate 
calculation of the radiation field. 

We avoid this inefficien cy by utilising adaptive ray trac- 
ing (|Abel Wandeltll2002l V which progressively splits rays 
when the sampling becomes too coarse and is based on 
Hierarchical Equa l Area isoLatitude Pixelation (HEALPix; 
iGorski et al.ll200^ ). In this scheme, the rays are traced along 
normal directions of the centres of HEALPix pixels, which 
evenly divides a sphere into equal areas. The rays are ini- 
tialised at each point source with the photon luminosity (ph 
s~^) equally spread across A^pix = 12 x 4^ rays, where / is 
the initial HEALPix level. We usually find / = or 1 is suf- 
ficient because these coarse rays will usually be split before 
traversing the first cell. 

The rays are traced through the grid in a typical fashion 
(e.g. lAbel et al.l[T999 ^. in which we calculate the next cell 
boundary crossing. The ray segment length crossing the cell 



dr — Ro 



min [(xceii,i 



c,i)/nr 



(5) 



where i?o, '^ray, a:ceii,i, and Xsrc,i are the initial distance trav- 
elled by the ray, normal direction of the ray, the next cell 
boundary crossing in the i-ih dimension, and the position of 
the point source that emitted the ray, respectively. However 
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Table 1. Variable definitions used in Sections [2] and [3] 

Variable Description 

Aceii Cell face area 

a Scale factor 

Ch Collisional ionization rate of hydrogen 

CRT,cfl CFL safety factor in timestep calculation 

Dc,i Distance from ray segment center to cell center 

Dedge Distance from ray segment center to cell edge 

dP Photon loss from absorption 

dPc Photon loss from Compton Scattering 

dp-y Momentum change from radiation pressure 

dtp Photon timestep 

£^ph Photon energy 

Ei Ionization energy of absorber 

fc Geometric correction factor 

/shield Shielding function for H2 

H Hubble constant 

Iv Specific intensity 

Emission coefficient 

A^ph Photo-ionization rate 

^diss Photo-dissociation rate of H2 

Lpix Linear width of a HEALPix pixel 

/ HEALPix level 

A^H2 Column density of H2 

''^abs Number density of absorber 

^pix(0 HEALPix pixels on level I 

Nray Rays per cell 

A^^ Photon luminosity 

n Normal direction of radiation 

Ux Number density of absorber x 

P Photon ffux 

r Radius 

Tceii Distance from the radiation source to the 
cell center 

Tray Distance from the radiation source to the 

next cell boundary crossing 

Vceii Cell volume 

viF Ionization front velocity 

xo,i Cell center coordinates 

3^cell,i Next cell boundary crossing in the 

i-th dimension 

Xsrc,i Source position in the i-th dimension 

Secondary ionization factors 

aB Case-B recombination rate 

Fph Photo-heating rate 

Ax Cell width 

<I>c Angular resolution in units of rays per cell area 

hiu Absorption coefficient 

^mfp Mean free path 

r^ray Solid angle associated with a ray 

(^ahs Cross-section of absorber 

^ray Angle associatcd with a ray 

r Optical depth 



before the ray travels across the cell, we evaluate the ratio 
of the face area Acew of the current cell and the solid angle 



r^ray of the ray, 



A,, 



iVpix(Aa;)' 



(6) 



If $c is less than a pre-determined value (usually > 3), the 
ray is split into 4 child rays. We investigate the variations 
in solutions with $c in ^7.21 The pixel numbers of the child 
rays p' are given by the "nested" scheme of HEALPix at 



the next level, i.e. ^ = 4 x p + [0, 1,2,3], where p is the 
original pixel number. The child rays (1) acquire the new 
normal vectors of the pixels, (2) retain the same radius of 
the parent ray, and (3) gets a quarter of the photon flux of 
the parent ray. Afterwards the parent ray is discontinued. 
A ray propagates and splits until 

(i) the photon has travelled c x dtp, where dtp is the 
radiative transfer timestep, 

(ii) its photon flux is almost fully absorbed (> 99.9%) in 
a single cell, which significantly reduces the computational 
time if the radiation volume filling fraction is small, 

(iii) the photon leaves the computational domain with 
isolated boundary conditions, or 

(iv) the photon travels of the simulation box length 
with periodic boundary conditions. 

In the first case, the photon is halted at that position and 
saved, where it will be considered in the solution of Ijy at 
the next timestep. In the next timestep, the photon will 
encounter a different hydrodynamical and ionisation state, 
hence in its path. Furthermore any time variations of the 
luminosities will be retained in the radiation field. This is 
how this method retains the time derivative of the radiative 
transfer equation. The last restriction prevents our method 
from considering sources external to the computational do- 
main, but a uniform radiation background can be used in 
conjunction with ray tracing in ENZO+MORAY that adds the 
local radiation field to the background intensity. 



2.2 Radiation Field 

The radiation field is calculated by integrating equation 
dH) along each ray, which is done by considering the dis- 
cretisation of the ray into segments. In the following sec- 
tion, we assume the rays are monochromatic. For conve- 
nience, we express the integration in terms of optical depth 
T = J i^(r,t) dr, and for a ray segment, 

dr = crabs(z^)nabsC?r. (7) 

Here dabs and ^abs are the cross section and number den- 
sity of the absorbing medium, respectively. We use the cell- 
centred density in our calculations. Using trilinearly inter- 
polated densities (see Meller oa et al.i r20Q6) did not produce 
improved results. In the static case, equation (|4| has a sim- 
ple exponential analytic solution, and the photon flux of a 
ray is reduced by 



dP = Px (1-e"^) 



(8) 



as it crosses a cell. We equate the photo-ionisation rate 
to the absorption rate, resulting in p hoton conservation 
(|Abel et al.1 Il999l : iMellema et "alTl2006l ). Thus the photo- 
ionisation kph and photo-heating Fph rates associated with 
a single ray are 



ph 



P(l 



riabs VceW dtp 

Fph = kph {Eph — Ei), 



(9) 
(10) 



where VceU is the cell volume, ^ph is the photon energy, 
and Ei is the ionisation energy of the absorbing material. 
In each cell, the photo-ionisation and photo-heating rates 
from each ray in the calculation are summed, and after the 
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Figure 1. (a) A two-dimensional illustration of the overlap be- 
tween the beam associated with a ray 7 and a computational 
cell. The ray has a segment length of dr passing through the cell. 
The covering area is denoted by dark grey, where the full area 
(dr X Lpix) is colored by the dark and light grey. The photo- 
ionisation and photo-heating rates should be corrected by the 
ratio of these areas, i.e. the overlap fraction fc. (b) Annotation 
of quantities used in this geometric correction. 



ray tracing is complete, these rates can be used to update 
the ionisation state and energy of the cells. Considering a 
system with only hydrogen photo-ionisations and radiative 
recombinations, these changes are very straightforward and 
is useful for illustrative purposes. The change in neutral hy- 
drogen is 

dnu 

~dr 

where as = 2.59 x 10~^^cm^ s~^ is the recombination co- 
efficient at 10^ K in the Case B on-the-spot approximation 
in w hich all recombinations are locally reabsorbed, (Spitzer 
[l978), and Ch is the collisional ionisation rate. However, for 
more accurate solutions in calculations that consider several 
chemical species, the photo-ionis ation rates are te rms in the 
relevant chemical networks fe.g. lAbel et al.|[l997i V 



aBfleflp — CHfleTlH — k 



(11) 



2.3 Geometric Corrections 

For a ray tracing method to accurately, i.e. without non- 
spherical artifacts, compute the radiation field, the compu- 
tational grid must be well-sampled by the rays. The main 
source of potential artifacts is the geometrical difference be- 
tween the cell and the HEALPix pixel relevant in the angular 
integration of the intensity over the cell. In this section, we 
devise a correction scheme to account for these differences. 
Consider the solid angle Qray and photon flux P associated 
with a single ray, and assume the flux is constant across 
r^ray. Thcrc cxists a discrepancy between the geometry cell 
face and HEALPix pixel when the pixel does not cover the 
entire cell face, which is illustrated in Fig. [1] This mismatch 
causes non-spherical artifacts and is most apparent in the 
optically thin case, where the area of the pixel is dominant 
over (1 — e~^) when calculating kph- One can avoid these ar- 
tifacts by increasing the sampling $c to high values (> 10) 
but we have formulated a simple geometric correction to 
the calculation of the radiation field. This correction is not 
unique to the HEALPix formalism but can be applied to 
any type of pixelisation. 

The contribution to kph and Fph must be corrected by 
the covering factor fc of the pixel with respect to the cell. 



When the pixel is fully contained within the cell face, fc = 
1. Because the geometry of the pixel can be complex with 
curved edges, we approximate fc by assuming the pixel is 
square. The covering factor is thus related to the width of a 
pixel, Lpix = Ro ^pix, and the distance from the ray segment 
midpoint to the closest cell boundary Dc,i, which is depicted 
in Fig. [1] To estimate /c, we first find the distance (icentre,i 
from the midpoint of the ray segment to the cell centre xo,i 
in orthogonal directions. 



RoA + rii 



dr 



(12) 



where i?o,i is the distance travelled by the ray in each or- 
thogonal direction. The distance to the closest cell boundary 
is Dedge = dx/2 — maxi=i^3(Dc,i)- Thus the covering factor 
is related to the square of the ratio between the Lpix and 
L^edge by 



fc 



1 

2+ L 



edge 



(13) 



One half of the pixel is always contained within the cell, 
which results in the factor of 1/2. Finally we multiply /cph 
and Fph by fc but leave the absorbed radiation dP un- 
touched because this would underestimate the attenuation of 
the incoming radiation. Using fc calculated like above, the 
method is no longer photon conserving. In our implemen- 
tation, we felt that the spherical symmetry obtained out- 
weighed the loss of photon conservation. However we show 
that there are no perceptible deviations from photon con- 
servation in ^4. II and ^7.11 

We briefly next describe how to retain photon conserva- 
tion with a geometric correction. Notice that we compute fc 
by only considering the distances in orthogonal directions. 
A better estimate would consider the distance between the 
cell boundary and ray segment midpoint in the direction 
between the midpoint and cell centre of Xmid — xq. We find 
that the method outlined here provides a sufficient correc- 
tion factor to avoid any non-spherical artifacts and devia- 
tions from photon conservation. Furthermore in principle, 
the ray should also contribute to any neighboring cells that 
overlap with Qray , which is the key to be photon conservative 
with such a geometric correction. 



2.4 Optically Thin Approximation 

In an optically thin medium, radiation is only attenuated by 
geometric dilution in the local approximation to Equation 
([2]), i.e. the inverse square law. With such a simple solution, 
the tracing of rays is wasteful, however these rays must be 
propagated because the the medium farther away can be 
optically thick. Here we describe a method that minimizes 
the computational work of ray tracing in the optically thin 
regime by exploiting this fact. Each ray tracks the total col- 
umn density A^abs and the equivalent total optical depth r 
traversed by the photon. If r < rthin ~ 0.1 after the ray 
exits the cell, we calculate the photo-ionisation and photo- 
heating rates directly from the incoming ray instead of the 
luminosity of the source. 



C^abs 



dtt 



^rav 



(14) 
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Note that the photon number P in the ray has already been 
geometrically diluted by ray splitting. Here rceii and Tray 
are the radii from the radiation source to the cell centre and 
where the ray exits the cell. Thus the last factor corrects the 
flux to a value appropriate for the cell centre. The photo- 
heating ray is calculated in the same manner as the general 
case, Fph = kph{Eph — Ei). This should only be evaluated 
once per cell per radiation source. No photons are removed 
from the ray. With this method, we only require one ray 
travel through each cell where the gas is optically thin, thus 
reducing the computational expense. 

We must be careful not the overestimate the radiation 
when multiple rays enter a single cell. In the case of a sin- 
gle radiation source, the solution is simple - only assign 
the cell the photo-ionisation and photo-heating rates when 
A:ph = 0. However in the case with multiple sources, this is 
no longer valid, and we must sum the flux from all optically 
thin sources. Only one ray per source must contribute to a 
single cell in this framework. We create a flagging field that 
marks whether a cell has already been touched by an opti- 
cally thin photon from a particular radiation source. Naively, 
we would be restricted to tracing rays from a single source 
at a time if we use a boolean flagging field. However we can 
trace rays for 32 sources at a time by using bitwise opera- 
tions on a 32-bit integer field. For example in C, we would 
check if an optically thin ray from the n-th source has prop- 
agated through cell i by evaluating (MarkerField[i] ^ n 
& 1). If false, then we can add the optically thin approxi- 
mation [equation (|14p ] to the cell and set MarkerField[i] 
1= (1 ^ n) ; to mark the cell. 

2.5 Additional Physics 

Other radiative processes can also be important in some 
situations, such as attenuation of radiation in the Lyman- 
Werner bands, secondary ionisations from X-ray radiation, 
Compton heating of from scattered photons, and radiation 
pressure. We describe our implementation of these physics 
next. 



2.5.1 Absorption of Lyman- Werner Radiation 

Molecular hydrogen can absorb photons in the Lyman- 
Werner bands through the two-step Solomon process, 
which for the lowest ro-vibrational states already con- 
sists of_76_absor2tion lines ranging from 11.1 to 13.6 
eV (IStecher Williamsl Il967l : Dakarno k Stephens! Il970l : 
iHaiman et al.1 l2000l ). Each of these spectral lines can be 
modelled with a ty pical exponential attenuation equation 
(jRicotti et al.ll200lh . but Draine k Bertoldi (1996) showed 
that this self-shielding is well modeled with the following 
relation to total H2 column density 



/shield (A^H2) = 



1 (iVH2 ^ 10^^) 

(iVH2/10^^)-°-^^ (iVH2 > 10^^) 



(15) 



where A^h2 is in units of cm~^. To incorporate this shielding 
function into the ray tracer, we store the total H2 column 
density and calculate the H2 dissociation rate by summing 
the contribution of all rays. 



kd 



P CTLW ^ray 

rays 



^cell dV dtp 



(16) 



whe re ctlw — 3 . 71 x IQ^^cm^ is the effective cross-section of 
H2 (jAbel et al.lll997l l. To account for absorption, we atten- 
uate the photon number flux by 

dP = P[/shield(iVH2 + dN^2) - /shield(iVH2)], (17) 

where c/A^h2 is the H2 column density in the current cell. 

2.5.2 Secondary Ionisations from X-rays 

At the other end of the spectrum, a high-energy (^ph ^ 100 
eV) photon can ionise multiple neutral hydrogen and he- 
lium atoms, an d this should be con sidered in such radiation 
fields. iShull van Steenberg (1985) studied this effect with 
Monte Carlo calculations over varying electron fractions and 
photon energies up to 3 keV. They find that the excitation 
of hydrogen and helium and the ionisation of He 11 is negli- 
gible. The number of secondary ionisations of H and He is 
reduced from the ratio of the photon and ionisation energies 
(Eph/Ei) by a factor of 



yk,H = 0.3908(1 



0.4092N 1.7592 
X ) 



yk,He = 0.0554(1 - X° 



(18) 
(19) 



where x is the electron fraction. The remainder of the photon 
energy is deposited into thermal energy that is approximated 
by 



Yr = 0.9971[1 - (1 



0.2663\1.3163i 



(20) 



and approaches one as x ^ 1. Thus in gas with low electron 
fractions, most of the energy results in ionisations of hydro- 
gen and helium, and in nearly ionised gas, the energy goes 
into photo-heating. 

2. 5. 3 Compton Heating from Photon Scattering 

High energy photons can also cause Compton heating by 
scattering off free electrons. During a scattering, a photon 
loses AE{Te) = 4A:Te x (Eph/meC^) of energy, where Te is 
the electron temperature. For the case of monochromatic 
energy groups, we model this process by considering that 
the photons are absorbed by a factor of 



dPc 



-r.^A(T; 



(21) 



which is the equivalent of the photon energy decreasing. Here 
Te = UeO-KNdl is the optical depth to Compton scattering, 
and (Jkn is the non-relativ istic Klein-Nishina cross section 
(jRvbicki Lightmanlll979h . The Compton heating rate is 
thus 



- ph,C — 



dPc 



(22) 



Tie Vcell dtp ' 

which has been used in lKim et al.l (|201ll ). 
2.5.4 Radiation Pressure 

Another relevant process is radiation pressure, where the 
absorption of radiation transfers momentum from photons 
to the absorbing medium. This is easily computed by con- 
sidering the momentum 



dpj 



dP ^ph 



(23) 
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of the absorbed radiation from the incoming ray, where f is 
the normal direction of the ray. We do not include radiation 
pressure on dust currently. The resulting acceleration of the 
gas because of radiation pressure is 

rfa= ./P?;, . (24) 

dtp p Vcell 

where p is the gas density inside the cell. This acceleration is 
then added to the other forces, e.g. gravity, in the calculation 
in an operator split fashion. 



3 NUMERICAL IMPLEMENTATION IN ENZO 

In this section, we describe our parallel implementation of 
the adaptive ray tracing m ethod into e nzo. enzo is a paral- 
lel block-structured A MR (iBerger Co lella 1989 ) code that 
is pu blicly available (iBrvan &: Norman 1997; O 'Shea et al.l 
First we explain the programming design of handling 
the "photon packages" that are traced along the adaptive 
rays. We use the terms photon packages and rays inter- 
changeably. Next we focus on the details of the radiation 
hydrodynamics and then the importance of correct time- 
stepping. Last we give our parallelisation strategy of tracing 
rays through an AMR hierarchy. This implementation is in- 
cluded in the v2 . public version of enzo. 

3.1 Programming Design 

Each photon package is stored in the AMR grid with the 
finest resolution that contains its current position. The pho- 
ton packages keep track of their (1) photon flux, (2) photon 
type, (3) photon energy, (4) the time interval of its emission, 
(5) emission time, (6) current time, (7) radius, (8) total col- 
umn density, (9) HEALPix pixel number, (10) HEALPix 
level, and (11) position of the originating source, totaling 
60 (88) bytes for single (double) precision. When enzo uses 
double precision for grid and particle positions and time, 
items 4-7 and 11 are double precision. 

We only treat point sources of radiation in our imple- 
mentation; therefore all base level photon packages origi- 
nate from them. As they travel away from the source, they 
generally pass through many AMR grids, especially if the 
simulation has a high dynamic range. This is a challenging 
programming task as rays are constantly entering and ex- 
iting grids. Before any computation, the number of rays in 
a particular grid is highly unpredictable because the inter- 
vening medium is unknown. Furthermore, the splitting of 
parent rays into child rays and a dynamic AMR hierarchy 
add to the complexity. Because of th is, we store th e pho- 
ton packages as a doubly hnked list ((Abel k Wandeltir2002l ) . 
Thus we can freely add and remove them from grids without 
the concern of allocating enough memory before the tracing 
commences. 

We illustrate the underlying algorithm of the ray tracing 
module in enzo in Fig. [2] and the ray tracing algorithm is 
shown in Fig. [S] The module is only called when advancing 
the finest AMR level. We describe its steps below. 

Step 1- Create A^pix new photon packages on the initial 
HEALPix level from point sources. Place the new rays in the 
highest resolution AMR grid that contains the source. 

Step 2- Initialise all radiation fields to zero. 



Step 3- Loop through all AMR grids, tracing any rays 
that exist in it. For each ray, the following substeps are 
taken. 

Step 3a- Compute the ray normal based on the 
HEALPix level and pixel number of the photon package 
with the HEALpix routine pix2vec_nest. One strategy to 
accelerate t he computation is to store ray segment path s 
in memory (|Abel &: Wandeltl I2OO2I : iKrumholz et al.ll2007l ): 
however this must be recomputed if the grid structure or 
point source position changes. We do not restrict these two 
aspects and cannot employ this acceleration method. 

Step 3b- Compute the position of the ray (rsrc + rn), 
the current cell coordinates in floating point and its cor- 
responding integer indices. Here rsrc is the position of the 
point source, r is the distance travelled by the ray, and n is 
the ray normal. 

Step 3c - Check if a subgrid exists under the current 
cell. If so, move the ray to a linked list that contains all 
rays that should be moved to other grids. We call this vari- 
able PhotonMoveList. Store the destination grid number and 
level. Continue to the next ray in the grid (step 3a). We de- 
termine whether a subgrid exists by creating a temporary 
3D field of pointers that either equals the pointer of the cur- 
rent grid if no subgrid exists under the current cell or the 
child grid pointer that exists under the current cell. This 
provides a significant speedup when compared to a simple 
comparison of a photon package position and all of the child 
grid boundaries. Note that this is the same algorithm used 
in ENZO when moving collisionless particles to child grids. 

Step 3d - Compute the next cell crossing of the ray and 
the ray segment length across the cell (Equation [5|). 

Step 3e- Compare the solid angle associated with the 
ray at radius r -\- dr with a user-defined splitting criterion 
(Equation [6]). If the solid angle is larger than the desired 
minimum sampling, split the ray into 4 child rays f ^2.1|) . 
These rays are inserted into the linked list after the parent 
ray, which is subsequently deleted. Continue to the next ray 
(step 3a), which will be the first child ray. 

Step 3f- Calculate the geometric correction (Equation 
[13]), the optical depth of the current cell (Equation[7|), photo- 
ionisation and photo- heating rates (Equations [9] and [TO)) . 
and add the column density of the cell to the total column 
density of the ray. 

Step 3g.- Add the effects of any optional physics mod- 
ules f ^2.5|) - secondary ionisations from X-rays, Compton 
heating from scattering, and radiation pressure. 

Step 3h.- Update the current time (t = t + cdr), photon 
flux {P — P — dP, Equation [8]) , and radius of the ray (r = 
r + dr). 

Step 3i- If the photon flux is zero or the total optical 
depth is large (> 20), delete the ray. 

Step 3j- Check if the ray has travelled a total distance 
of cdtp in the last timestep. If we are keeping the time- 
derivative of the radiative transfer equation, halt the photon. 
If not (i.e. infinite speed of light), delete the photon. 

Step 3k.- Check if the ray has exited the current grid. 
If false, continue to the next cell (step 3b). If true, move 
the ray to the linked list PhotonMoveList, similar to step 
3c. If the ray exits the simulation domain, delete it if the 
boundary conditions are isolated; otherwise, we change the 
source position of the ray by a distance -sign(n[i] ) * 
DomainWidthEi] , where n is the ray normal, and i is the 
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Create base photons 
from pt. sources 




EXIT to 
main grid loop 



Figure 2. Flow chart for the overall algorithm of the radiative transfer module in ENZO that illustrates (1) the creation of photon 
packages, (2) ray tracing, (3) the transport of photon packages between AMR grids, and (4) coupling with the hydrodynamics. The ray 
tracing algorithm, which is contained in the "Trace Rays" is detailed in Fig. [S] 



dimension of the outer boundary it has crossed. The radius 
is kept unchanged. In essence, this creates a "virtual source" 
outside the box because the ray will be moved to the oppo- 
site side of the domain, appearing that it originated from 
this virtual source. 

Step 4~ If ciny rays exist in the linked list 
PhotonMoveList, move them to their destination grids and 
return to step 3. This requires MPI communication if the 
destination grid exists on another processor. 

Step 5- If all rays have not been halted (keeping the 
time- derivative of the radiative transfer equation), absorbed, 
or exited the domain, return to step 3. 

Step 6- With the radiation fields updated, call the 
chemistry and energy solver and update only the cells with 
radiation, which is discussed further in ^3.31 

Step 7- Advance the time associated with the photons 
tp by the global timestep dtp (for its calculation, see ^3.4|) . If 
tp does not exceed the time on the finest AMR level, return 
to step 1. 



3.2 Energy groups 

In our implementation, p hoton packages are m ono- 
chromatic, i.e. energy groups (|Mihalas Mihalaslll984l . Ch. 
6) , and are assigned a photon type that corresponds whether 
it is a photon that (1) ionises hydrogen, (2) singly ionises 
helium, (3) doubly ionises helium, (4) has an X-ray energy, 
or (5) dissociates molecular hydrogen (Lyman- Werner ra- 
diation). One disadvantage of mono-chromatic rays is the 
number of rays increase with the number of frequency bins. 
However this allows for early termination of rays that are 
fully absorbed, which are likely to have high absorption 
cross-sections (e.g. H i ionisations near 13.6 eV) or a low 



initial intensity (e.g. He ii ionising photons in typical stel- 
lar populations) . The other approach used by some groups 
fe.g. lTrac fc"c^l2Q07l ) is to store all energy groups in a sin- 
gle ray. This reduces the number of the rays generated and 
the computation associated with the ray tracing. Unless the 
ray dynamically adjusts its memory allocation for the en- 
ergy groups as they become depleted, this method is also 
memory intensive in the situation where most of the energy 
groups are completely absorbed but a few groups still have 
significant flux. 

In practice, we have found that one energy group per 
photon type is sufficient to match expected analytical tests 
f ^7.3|) . For exampl e when modeling Population HI stellar 
radiation (e.g. A bel et ahllioOTl : IWise k A bel 2008b, for hy- 
drogen ionising radiation only), we have 3 energy groups - 
H I, He I, He ll - each with an energy that equals the average 
photon energy above the ionisation threshold. 

3.3 Coupling with Hydrodynamics 

Solving the radiative transfer equation is already an inten- 
sive task, but coupling the effects of radiation to the gas 
dynamics is even more difficult because the radiation fields 
must be updated on a time-scale such that it can react to the 
radiative heating, i.e. sound-crossing time. The frequency of 
its evaluation will be discussed in the next section. 

ENZO solves the physical equations in an operator-split 
fashion over a loop of AMR grids. On the finest AMR level, 
we call our radiation transport solver before this main grid 
loop in the following sequence: 

• All grids: 

(i) Solve for the radiation field with the adaptive ray 
tracer 
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Pre-compute ray normal, ray position, 
current cell, and cross-section 



Does a child grid 
exist under this cell?. 



YES ^ 


Put photon in move list 




EXIT 



Compute next cell crossing 



YES 



Does the ray need 
splitting at r+dr? 



NO 




Split ray. Delete parent ray 
EXIT 



Calculate 
1 . Geometric correction 
2. Optical depth 
3. Photo-ionisation and photo-heating rates 
4. Add to column density 



Update cell position 



NO 



Optional: add radiation pressure 




If constant timestep, halt photon. 
If infinite c, delete. 
EXIT 



NO 



Update photon time, flux, and radius 




. YES^ 


Delete photon 




EXIT 



Figure 3. Flow chart for the ray tracing algorithm for one photon passing through a grid. Note that only one step is needed in the 
routine to adaptively split rays. The remainder is a typical ray tracing method. 



(ii) Update species fractions and energies for cells with 
radiation with a non-equilibrium chemistry solver on sub- 
cycles (Equation [25]). 

• For each grid: 

(i) Solve for the gravitational potential with the particle 
mesh method 

(ii) Solve hydrodynamics 

(iii) Update species fractions and energies for cells with- 
out radiation with a non-equilibrium chemistry solver on 
subcycles (Equation [25]). 

(iv) Update particle positions 

(v) Star particle formation 

• All grids: Update solution from children grids 

Since the solver must be called many times, the effi- 
ciency of the radiation solver is paramount. After every ra- 
diation timestep, we call the non-equilibrium chemistry and 
energy solver in enzo. This solves both the energy equa- 
tion and the network o f stiff chemical equa tions on small 
timesteps, i.e. subcycles (jAnninos et al.|[l997l V The timestep 
is 

at - mm |^^^^^^| , ' ' 2 ; ' ^ ^ 



where Ue is the electron number density, e is the specific 
energy, and hydro is the hydrodynamic timestep. This lim- 
its the subcycle timestep to a 10% change in either elec- 
tron density, neutral hydrogen density, or specific energy. 
In simulations without radiation, enzo calls this solver in a 
operation-split manner after the hydrodynamics module for 
grids only on the AMR level that is being solved. In simula- 
tions with radiative transfer, the radiation field can change 
on much faster time-scales than the normal hydrodynamical 
timesteps. 

For example, a grid on level L might have no radiation 
in its initial evaluation, but the ionisation front exists just 
outside its boundary. Then radiation permeates the grid in 
the time between tL ti,-\-dti,^ and the energy and chemical 
state of the cells must be updated with each radiation up- 
date to advance the ionisation front accurately. If one does 
not update these cells, it will appear that the ionisation 
front does not enter the grid until the next hydrodynami- 
cal timestep! Visually this appears as discontinuities in the 
temperature and electron fraction on grid boundaries. One 
may avert this problem by solving the chemistry and energy 
equations for every cell on every radiative transfer timestep, 
but this is very time consuming and unnecessary, especially 
if the radiation filling factor is small. 
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We choose to dynamically split the problem by cells 
with and without radiation. In every radiation timestep, the 
chemo-thermal state of only the cells with radiation are up- 
dated. For the solver subcycling, we replace c/thydro with dtp 
in Equation [25] in this case. Once the radiative transfer solver 
is finished with its timesteps, the hydrodynamic module is 
called, and then the chemo-thermal state of the cells with- 
out radiation are updated on a subcycle timestep stated in 
Equation [25] 

For cells that transition from zero to non-zero photo- 
ionisation rates, the initial state that enters into the chem- 
istry and energy solver does not correspond to the current 
time of the radiation transport solver tuT , but either time tL 
if the grid level is the finest level because its chemo-thermal 
state has not been updated or time tL + dti, on all other lev- 
els. In principle, one could first revert the cell back to time 
tL and then update to tRx with the chemistry and energy 
solver if the cell is on the finest level. However in practice, 
the time-scales in gas without radiation are small compared 
to the ionisation and heating time-scales when radiation is 
introduced. Therefore, we do not perform this correction and 
find that this does not introduce any inaccuracies in both 
test problems (see 311) and real world applications. 



3.4 Temporal evolution 

There have been several methods of choosing a maximal 
timestep to solve radiation transfer equation while retaining 
stability and accuracy. We describe several methods to cal- 
culate the radiative transfer timestep in this section. With a 
small enough timestep, the solution is accurate (ignoring any 
systematic ones), but the solver is slow. Furthermore for very 
small timesteps, the photon packages only advance a short 
distance, and they will exist in every dx/dtp cells with radi- 
ation and are stored between timesteps, excessively consum- 
ing memory. On the shortest time-scale, one ca n safely set 
the t i mestep to the ligh t- crossing time of a cell (lAbel et al.l 
ll999l : lTrac Cenll2Q07l l but encounters the problems stated 
above. 

If the timestep is too large, the solution will become 
inaccurate; specifically, ionisation fronts will advance too 
slowly, as radiation intensity exponentially drops with a 
scale length of the mean free path 



A 



1 



mfp = (26) 

past the ionisation front. For example in our implementa- 
tion, the chemo-thermal state of the system remains con- 
stant as the rays are traced through the cells. In the case 
of a single H ii region, the speed of the ionisation front is 
limited to approximately Xmip/dtp. 



3.4' i Minimizing neutral fraction change 

Another strategy is restricting the neutral fraction to change 
a small amount, i.e. for a single cell. 



dt 



nm 



^ 'il^H ^ (27) 

\dnm/dt\ |/cph + ne(CH + as)| ' 

wher e eion is the maximu m fraction change in neutral frac- 
tion. I Shapiro et al.l ()2004l ) found that this limited the speed 
of the ionisation front. We can investigate this further 




J 0.0020 




Figure 4. Radiative transfer adaptive timestep in shadowing test 
(Test 3; ^4.3|) while restricting the neutral fraction change to 5% 
in the ionisation front. The unmodified timestep (top) is slightly 
more noisy and the minima are more prominent than the timestep 
computed with a running average of the last two timesteps (bot- 
tom). The points show every tenth timestep taken into account 
for the running average. The sawtooth behavior is created by the 
ionisation front advancing into the next neutral cell in the over- 
density. 



by evaluating the ionisation front velocity in a growing 
Stromgren sphere without recombinations, 

= ATiR'^nHViF. (28) 

Using kph — N^a /A-kR^ AceW and /cph oc nn/nn, we can make 
substitutions respectively on both sides of Equation (|28)) to 
arrive at the ionisation front velocity viy oc nn/riYi. 

We have implemented this method but we only consider 
cells within the ionisation front (by experiment we choose 
r > 0.5) because we are interested in evolving ionisation 
fronts at the correct speed. In the ionised region, the abso- 
lute changes in neutral fraction are small and will not signif- 
icantly affect the ionisation front evolution. In other words, 
rtni / [duYii / dt) may be large but (dnm/dt) ^ 0, thus we 
can safely ignore these cells when determining the timestep 
without sacrificing accuracy. 

We search for the cell with the smallest dtp based on 
Equation 1271 In principle, one could use this value without 
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modifications as the timestep, but there is considerable noise 
both spatially and temporally. In order to stabilise this tech- 
nique, we first spatially smooth the values of dtp, cell with a 
Gaussian filter over a 3^ cube. Because we only consider the 
cells within the ionisation front, we set dtp^ceii to the hydro- 
dynamical timestep outside the front during the smooth- 
ing. After we have smoothed dtp, cell, we select the minimum 
value as dtp. Significant noise in dtp can exist between time 
steps as well. Because the solution can become inaccurate if 
the timestep is allowed to be too large, we restrict the next 
timestep to be less than twice the previous timestep, 

(itp,i = min(2(itp,o, dtp,i). (29) 

Fig. m shows the smooth evolution of dtp in a grow- 
ing Stromgren sphere when compared to the values of 
min((itp,ceii). 

3.4-2 Time averaged quantities within a timestep 

iMellema et al.l ([2006 ) devised an iterative scheme that al- 
lows for large timesteps while retaining accuracy by consid- 
ering the time- averaged values (r, /cph, rie, nm) during the 
timestep. Starting with the cells closest to the source, they 
first calculate the column density to the cell. Then they com- 
pute the time-averaged neutral density for the cell and its 
associated optical depth, which is added to the total time- 
averaged optical depth. With these quantities, one can com- 
pute a photo-ionisation rate and update the electron density. 
This process is repeated until convergence is found in the 
neutral number density. In a test with a Stromgren sphere, 
they found analytical agreement with 10""^ less timesteps 
than a method without time- averaging. Another advantage 
of this method is the use of pre-calculated tables for the 
photo-ionisation rates as a function of optical depth, based 
on a given spectrum. This minimizes the energy groups 
needed to accurately sample a spectrum. We are currently 
implementing this method into ENZO+MORAY. 



I199Q| \ the I- front detaches from the shock front, accelerates, 
and transitions back to an R-type front. This can also occur 
in champagne flows when the ionisation front passes a den- 
sity discontinuity. The I-front velocities in these two stages 
differ up to a factor of ^ 10. Although the solution is accu- 
rate with a large timestep in the D-type phase, the I-front 
may lag behind because of the constant timestep. After a few 
recombination times, the numerical solution eventually ap- 
proaches the analytical solution. If such a simulation focuses 
on the density core expansion and any small scale structures, 
such as cometary structures and photo-dissociation regions, 
one can cautiously sacrifice temporal accuracy at large scales 
for computational savings. 



3.4-4 Change of incident radiation 

Ionisation front velocities can approach significant fractions 
of the speed of light in steep density gradients and in the 
early expansion of the H ii region. If the ionisation front 
position is critical to the calculation, the radiation transport 
timestep can be derived from a non-relativistic estimate of 
the ionisation front velocity 



viF{r) 



F{r) ^ 
nabs(r) ' 



(30) 



based on the incident radiation field at a particular posi- 
tioiE. Alternatively, the propagation of the ionisation front 
can be restricted by limiting the change in specific intensity 
/ to a safety factor CRT,cfl, resulting in a timestep of 



dtp — Crta 



\dl/dt\' 



(31) 



We consider the specific intensity after the ray travels 
through the cell, so / = /oexp(— r), where r = nu(jdl is 
the optical depth through the cell. The time derivative of / 



3.4-3 Physically motivated 

A constant timestep is necessary when solving the time- 
dependent radiative transfer equation in enzo+moray. It 
should be small enough to evolve ionisation fronts accu- 
rately, as discussed earlier. The timestep can be based on 
physical arguments, for example, the sound-crossing time 
of an ionised region at T > 10^ K. To be conser vative, one 
may c hoose the sound-cros sing time of a cell fe.g. lAbel et al. 
I2OO7I : IWise Abell[2008bl ). Alternatively, the diameter of 
the smallest relevant system (e.g., an accretion radius, tran- 
sition radius to a D-type ionisation front, etc.) in the simu- 
lation may be chosen to calculate the sound-crossing time. 

If the timestep is too large, the ionisation front will 
propagate too slowly, but it eventually approaches the cor- 
rect radius at late times (see ^7.4p . This does not prevent one 
from using a large timestep, particularly if the system is not 
critically affected by a slower I-front velocity. One example is 
an expanding H 11 region in a power- law density gradient. Af- 
ter a brief, initial R-type phase, the I-front becomes D-type 
phase, where the ionisation and shock front progress jointly 
at the sound speed of the ionised region. A moderately large 
(0.1 Myr) timestep can accurately follow its ev olution. How- 
ever after the I-front passes a critical radius ([Franco et al.l 



^ — /o exp(-T)(-nHcrG 



(32) 



which can be expressed in terms of local optical depth and 
neutral fraction. 



dJ 
dt 



(33) 



flu is computed with the same formula as equation (|27)) . 
Substituting in equation (|3T)) gives 



dtp — CRT,cfl 



(34) 



In practice, we have found that a ceiling of 3 can be placed 
on the optical depth, so optically thick cells do not create a 
very small timestep. We still find excellent agreement with 
analytical solutions with this approximation. We show the 
accuracy using this timestep method in ^7.41 



^ See IShapiro et al ] ('2006') for the exact calculation of a rela- 
tivistic ionisation front. Neglecting relativistic terms do not affect 
the solution because the front velocity is only considered in the 
timestep calculation. 
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3.5 Parallelisation Strategy 

Parallelisation of the ray tracing code is essential when ex- 
ploring problems that require high resolution and thus large 
memory requirements. Furthermore, ENZO is already paral- 
lelised and scalable to O(IO^) processors in AMR simula- 
tions, and O(IO^) in unigrid calculations. ENZO stores the 
AMR grid structure on every processor, but only one pro- 
cessor contains the actual grid and particle data and photon 
packages. All other processors contain an empty grid con- 
tainer. As discussed in Step 4 in ^3.1[ we store the photon 
packages that need to be transferred to other grids in the 
linked list PhotonMoveList. In a single processor (serial) run, 
moving the rays is trivial by inserting these photons into the 
linked list of the destination grid. For multi-processor runs, 
we must send these photons through MPI communication 
to the processors that host the data of the destination grids. 
We describe our strategy below. 

The easiest case is when the destination grid exists on 
the same processor as the source grid, where we move the 
ray as in the serial case. For all other rays, we organise the 
rays by destination processors and send them in groups. We 
also send the destination grid level and ID number along 
with the ray information that is listed at the beginning of 

EH 

For maximum overlap of communication and computa- 
tion, which enables scaling to large numbers of processors, 
we must employ "non-blocking" MPI communication, where 
each processor does not wait for synchronisation with other 
processors. We use this technique for the sending and re- 
ceiving of rays. Here we desire to minimize the idle time of 
each processor when it is waiting to receive data. In the loop 
shown in Fig. [2] with the conditional that checks whether we 
have traced all of the rays, we aggressively transport rays 
that are local on the processor, and process any MPI re- 
ceive calls as they arrive, not waiting for their completion 
in order to continue to the next iteration. We describe the 
steps in this algorithm next. 

Step 1- Before any communication occurs, we count 
the number of rays that will be sent to each processor. The 
MPI receive calls (MPI.Irecv) must have a data buffer that 
is greater than or equal to the size of the message. We choose 
to send a maximum of A^max (= 10^ in ENZO v2 . 0) rays per 
MPI message. Therefore, we allocate a buffer of this size 
for each MPI.Irecv call. We then determine the number of 
MPI messages A^mesg and send this number in a non-blocking 
fashion, i.e. MPI_Isend. 

Step 2- Pack the photon packages into a contiguous 
array for MPI communication while the messages from Step 
1 completes. 

Step 3- Process the number of photon messages that 
we are expecting from each processor, sent in Step 1. Then 
post this number of MPI.Irecv calls for the photon data. Be- 
cause we strive to make the ray tracing routine to be totally 
non-blocking, the processors will most likely not be synchro- 
nised on the same loop (Steps 3-5 in ^3.1|) . Therefore, there 
might be additional Nmesg MPI messages waiting to be pro- 
cessed. We check for these messages and aggressively drain 
the message stack to determine the total number of photon 
messages that we are expecting and post their associated 
MPI.Irecv calls for the photon data. 



Step 4~ Send the grouped photon data with MPI.Isend 
with a maximum size of A^max photons. 

Step 5- Place any received photon data into the des- 
tination grids. We monitor whether the processor has any 
rays that were moved to grids on the same processor. If so, 
this processor has rays to transport, and we do not neces- 
sarily have to wait for any MPI receive messages and thus 
use MP I _Te St some to receive any messages that have already 
arrived. If not, we call MPI.Waitsome to wait for any MPI 
receive messages. 

Step 6- If all processors have exhausted their workload, 
then all rays have been either absorbed, exited the domain, 
or halted after travelling a distance cdtp. We check this in 
a similar non-blocking manner as the A^mesg calls in Step 1. 

Lastly we have experimented with a hybrid 
OpenMP/MPI version of ENZO, where workload is 
partitioned over grids on each MPI process. We found that 
parallelisation over grids for the photon transport does not 
scale well, and threading over the rays in each grid is a 
better approach. Because the rays are stored in a linked list 
in each grid, we must manually split the list into separate 
lists and let each thread work on each list. 



4 RADIATIVE TRANSFER TESTS 

Tests plays an important role in creating and maintaining 
computational tools. In this section, we present tests drawn 
from the Cosmological Radiative Transfer Codes Compar- 
ison Project (|lliev et al.l 1 20061 ), where results from 11 dif- 
ferent radiative transfer codes compared results in four 
test problems. The codes use various methods for radiation 
transport: ray tracing with short, long, and hybrid char- 
acteristics, Monte C arlo casting; ionisation front tracking 
( Alvarez et al.ll2006bll: variable Eddington Tensor formalism 
( Gnedin AbellboOlh . They conducted tests that investi- 
gated (1) the growth of a single Stromgren sphere enforcing 
isothermal conditions, (2) the same test with an evolving 
temperature field, (3) shadowing created by a dense, opti- 
cally thick clump, and (4) multiple H ii regions in a cos- 
mological density field. In all of the tests presented here, 
we use the method of restricted neutral fraction changes 
f ^3.4.ip for choosing a radiative transfer timestep. We cast 
48 rays (HEALPix level 1) from the point source and require 
a sampling of at least $c = 5.1 rays per cell. 



4.1 Test 1. Pure hydrogen isothermal H ii region 
expansion 

The expansion of an ionising region with a central source 
in a uniform m edium is a classic problem first studied by 
IStromgrenI ([l939i V This simple but useful test can uncover 
any asymmetries or artifacts that may arise from deficiencies 
in the method or newly introduced bugs in the development 
process. In this problem, the ionised region grows until re- 
combinations balance photo-ionisations [equation ([1])]. The 
evolution of the radius Vs and velocity Vs of the ionisation 
front has an exact solution of 



rs(t) = i?4l -exp(-t/tr, 



11/3 



Vs(t) = 



Rs 



exp(— t/tri 



3trec [1 -exp(-t/trec)]2/3' 



(35) 

(36) 
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Figure 5. Test 1 (H ll region expansion with a monochromatic 
spectrum of 13.6 eV). Top: Radially averaged profile of neutral 
(solid) and ionised (dashed) fraction at 10, 30, 100, and 500 Myr. 
Bottom: Evolution of the calculated (top, dashed) and analytical 
(top, solid) Stromgren radius. The ratio of these radii are plotted 
in the bottom panel. 



where tree = {aBnH)~^ is the recombination time. 

We adopt the problem parameters used in RT06. The 
ionising source emits 5 x 10^^ ph s~^ of monochromatic ra- 
diation at 13.6 eV and is located at the origin in a simula- 
tion box of 6.6 kpc. The ambient medium is initially set at 
T = 10"^ K, UH = 10"^ cm"^, X = 1.2 X 10"^, resulting in 
= 5.4 kpc and tree = 122.4 Myr. The problem is run for 
500 Myr. In the original tests, the temperature is fixed at 
10^ K; however, our solver is inherently tied to the chem- 
istry and energy solver. To mimic an isothermal behavior, 
we set the adiabatic index 7 = 1.0001, which ensures an 
isothermal state but not a fixed ionisation fraction outside 
of the Stromgren sphere. 

In Fig. O we show (a) the evolution of the neutral and 
ionisation fraction as a function of radius at t = 10, 30, 100, 
and 500 Myr, and (b) the growth of the ionisation front ra- 
dius as a function of time and its ratio with the analytical 
Stromgren radius [equation (|35|) ]. The ionisation front has a 
width of ^ 0.7 kpc, which is in agreement with the inherent 
thickness of ^ ISAmfp = 0.74 kpc, given a 13.6 eV mono- 
chromatic spectrum. There are small kinks in the neutral 
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Figure 6. Test 1 (H 11 region expansion with a monochromatic 
spectrum of 13.6 eV). Probability distribution function for neu- 
tral fraction at 10 Myr (solid), 100 Myr (dashed), and 500 Myr 
(dotted). Recombination of hydrogen at T = 10^ K causes the 
maximum neutral fraction to decrease from 1 to 0.75. 




Figure 7. Test 1 (H 11 region expansion with a monochromatic 
spectrum of 13.6 eV). Slice of neutral fraction at the origin at 500 
Myr. 



fraction at 1.5 and 3 kpc that corresponds to artifacts cre- 
ated by the photon package splitting at these radii. However 
these do not affect the overall solution. One difference be- 
tween our results and the codes presented in RT06 is the 
increasing neutral fraction outside of the H li region. This 
occurs because the initial ionised fraction and temperature 
is set tol.2xl0~^ and 10^ K, which are not the equilibrium 
values. Over the 500 Myr in the calculation, the neutral frac- 
tion increases to 0.2, which is close to its equilibrium value. 
In the right panel of Fig. \5\ the ionisation front radius ex- 
ceeds Rs by a few percent for most of the calculation. This 
difference happens because the analytical solution [equation 
(|35p ] assumes the H 11 region has a constant ionised frac- 
tion. The evolution of the ionised fraction as a function of 
radius can be analytically calculated (e.g. iOsterbrocklll989l : 



radms can be analytically 
IPetkova k Springell [2OO9I ) 



causing the ionisation front ra- 
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Figure 8. Test 2. (H ll region expansion with a T = 10^ K 
blackbody spectrum). Radially averaged profile of neutral (solid) 
and ionised (dashed) fraction at 10, 30, 100, and 500 Myr. 



10" 




-2.5 -2.0 -1.5 -1.0 -0.5 0.0 
log;LO Ionized Fraction 



2.0 2.5 3.0 3.5 4.0 4.5 

log^o Temperature (K) 



Figure 10. Test 2. (H ll region expansion with a T = 10^ K 
blackbody spectrum). Probability distribution function for ion- 
ized fraction (left) and temperature (right) at 10 Myr (solid), 100 
Myr (dashed), and 500 Myr (dotted). 
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Figure 9. Test 2. (H ll region expansion with a T = 10^ K 
blackbody spectrum). Evolution of the average neutral fraction. 



dius to be slightly larger, increasing from to 3% in the 
interval 80-350 Myr. Our results are in excellent agreement 
with this more accurate analytical solution. We show the 
distribution of neutral fraction in the domain for t — 10, 
100, and 500 Myr in Fig.[6]that agrees well with the results 
in RT06. In Fig. [71 we show a slice of the neutral fraction 
through the origin. Other than the ray splitting artifacts 
that generate the plateaus at 1.5 and 3 kpc, one sees spher- 
ical symmetry without any noise in our solution. 



4.2 Test 2. H ii region expansion: temperature 
evolution 

This test is similar to Test 1, but the temperature is al- 
lowed to evolve. The radiation source now has a blackbody 
spectrum with a T = 10^ K. The initial temperature is set 
at 100 K. The higher energy photons have a longer mean 
free path than the photons at the ionisation threshold in 
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Figure 12. Test 2. (H ll region expansion with a T = 10^ K black- 
body spectrum). Top: Evolution of the radius of the simulated 
ionisation front (dashed) and analytical (solid) Stromgren radius. 
Bottom: The ratio of the calculated and analytical Stromgren 
radius. 



Test 1. Thus the ionisation front is thicker as the pho- 
tons can penetrate deeper into the neutral medium. Here 
we use 4 energy groups with the following mean energies 
and relative luminosities: Ei = (16.74,24.65,34.49,52.06), 
Li/L = (0.277, 0.335, 0.2, 0.188). 

In Fig. [HI we show the radially averaged neutral and 
ionised fraction at t — 10, 30, 100, and 500 Myr, and the 
total neutral fraction of the domain in Fig. [9] Compared 
with Test 1, the ionisation front is thicker, as expected with 
the harder spectrum. Artifacts originating from ray split- 
ting, similar to Test 1, appear at r ~ 1 and 3 kpc as kinks 
in the neutral fraction. The total neutral fraction decreases 
to 0.67 over 4trec = 500 Myr, which is in agreement with 
the analytical expectation and other codes in RT06. Fig. [TOl 
shows the distribution of ionized fraction and temperature 
in Test 2. The overall trends agree with the codes presented 
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Figure 11. Test 2. (H ll region expansion with a T = 10^ K blackbody spectrum). Top: Slices through the origin of neutral fraction at 
10 and 100 Myr. Bottom: Slices of temperature at 10 and 100 Myr. 



in RT06 with the exception of the ray splitting artifacts that 
appear as slight rises in the distribution at logXe — 1 and 
logT ^ 3. In Fig. [TTl we show shces of neutral fraction and 
temperature through the origin at t = 10 and 100 Myr. Here 
one sees the spherically symmetric H ii regions and a smooth 
temperature transition to the neutral ambient medium. In 
Fig. 1121 we show the ratio of the ionisation front radius riF 
in our simulation and Rg. Before 1.5trec, Hf lags behind Rs, 
initially by 10% and then increases to Rs; however after- 
wards, this ratio asymptotes to a solution that is 4% greater 
than Rs. This behavior is approximately the median result 
in RT06, where this ratio varies between 1 and 1.1, and the 
early evolution of nr is under-predicted by almost all of the 
codes. If we use one energy group with the mean energy 
(29.6 eV) of a T = 10^ K blackbody, we find that rw/Rs = 
1.08, which is representative of the codes in the upper range 
of RT06. 



4.3 Test 3. I-front trapping in a dense clump and 
the formation of a shadow 

The diffusivity and angular resolution of a radiative trans- 
port method can be tested with the trapping of an ionisation 
front by a dense, neutral clump. In this situation, the ion- 
isation front will uniformly propagate until it reaches the 
clump surface. Then the radiation in the line of sight of the 
clump will be absorbed more than the ambient medium. If 
the clump is optically thick, a shadow will form behind the 
clump. The sharpness of the ionisation front at the shadow 
surface can be used to determine the diffusivity of the 
method. Furthermore the shadow surface should be aligned 
with the outermost neutral regions of the clump, which can 
visually assess the angular resolution of the method. 

The problem for this test is contained in a 6.6 kpc box 
with an ambient medium of uh — 2 x lO"'^ cm"*^ and T — 
8000 K. The clump is in pressure equilibrium with uh — 
0.04 cm~^ and T = 40 K. It has a radius of Tc = 0.8 kpc and 
is centred at (x, z) = (5, 3.3, 3.3) kpc. The ionized fraction 
of the entire domain is zero. In RT06, the test considered a 
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plane parallel radiation field with a flux Fq = 10^ ph 
cm~^ originating from the y = plane. Our code can only 
consider point sources, so we use a single radiation source 
located in the centre of the y — boundary. The luminosity 
of iV^ = 3 X 10^^ ph s~^ corresponds to the same flux Fq at 5 
kpc. The location wher e the ionisation front trapping can be 
calculated analytically ([Shapiro et al.ll2004l ). and with these 
parameters, it should halt at approximately the centre of 
the clump. We use the same four energy groups as in Test 
2. 

In Fig. 1131 we show neutral fraction and temperature 
of a one-dimensional cut through the centre of the dense 
clump at ^ = 0.5 at t = 1,3,5,15 Myr. The ionisation 
front is halted at a little more than halfway through the 
clump, which is consistent with the analytical expectation. 
The hardness of the T — 10^ K blackbody spectrum allows 
the gas outside of the ionisation front to be photo- heated. 
Where the gas is ionised, the temperature is between 10,000 
and 20,000 K, but decreases sharply with the ionised frac- 
tion. Fig. [13] depicts the average ionised fraction and tem- 



perature inside the dense clump, which both gradually in- 
crease as the ionisation front propagates through the over- 
density. Our results are consistent with RT06. The distribu- 
tion of ionized fraction and temperature within the clump 
are shown in Fig. [15] and agree well with the other codes in 
RT06. Finally we show slices of neutral fraction and temper- 
ature in the z = 0.5 plane in Fig. [TBI The neutral fraction 
slices prominently show the sharp shadows created by the 
clump and demonstrates the non-diffusivity behavior of ray 
tracing. The discretisation of the sphere creates one neutral 
cell on the left side of the sphere. This inherent artifact to 
the initial setup carries through the calculation. We did not 
smooth the clump surface like in some of the RT06 codes, in 
order to remove this artifact. It is seen in the neutral frac- 
tion and temperature states at all times and is not a caused 
by our ray tracing algorithm. 
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Figure 19. Test 4. (Multiple cosmological sources). Top: Slices through the origin of neutral fraction at 50 and 200 kyr at the coordinate 
z = 2;box/2. Bottom: Slices of temperature at 50 and 200 kyr. No smoothing has been applied to the images. 
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Figure 13. Test 3. (I- front trapping in a dense clump and shad- 
owing). Line cut from the point source through the middle of the 
dense clump at t = 1,3, 5, 15 Myr. of the average neutral fraction 
(top) and temperature (bottom) of the clump. 
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Figure 14. Test 3. (I-front trapping in a dense clump and shad- 
owing). Evolution of the average ionised fraction (top) and tem- 
perature (bottom) of the overdense clump. 
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Figure 15. Test 3. (I- front trapping in a dense clump and shad- 
owing). Probability distribution function for ionized fraction (left) 
and temperature (right) at 1 Myr (solid), 3 Myr (dashed), and 15 
Myr (dotted) within the dense clump. 
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Figure 17. Test 4. (Multiple cosmological sources). Probability 
distribution function for neutral fraction (left) and temperature 
(right) at 50 kyr (solid) and 200 kyr (dashed). 
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Figure 18. Test 4. (Multiple cosmological sources). Evolution of 
the mass- (dashed) and volume-averaged (solid) ionised fraction. 



agreement with the codes presented in RT06. We show 
the growth of the H ii regions by computing the mass- 
averaged Xm and volume- averaged Xv ionised fraction in 
Fig. 1181 Initially Xm is larger than Xv, and at t ^ 170 
kyr, the Xy beco mes larger. This is indicative of inside-out 
reionisation fe^. [ Gnedinl l2000l : iMiralda-Escude et al.1 l2000l : 
'Sokasian e t al.ll2003l v where the dense regions around haloes 
are ionised first, then the voids are ionised last. At the end 
of the simulation, 65% of the simulation is ionised. However 
by visual inspection in the slices of electron fraction (Fig. 
[19]), there appears to be very good agreement with C^-ray 
and FTTE. By first glance, our result appears to be different 
than the RT06 because of the color-mapping. Our results are 
also in good agreement with the ni ulti- frequency version of 
TRAPHIC (jPawlik Schavd 12011 see also for better rep- 
resentations of the electron fraction slices). In the slices of 
electron fraction and temperature. Fig. [191 the photo-heated 
regions are larger than the ionised regions by a factor of 2-3 
because of the hardness of the T ■ 
trum. 



10^ K blackbody spec- 



4.4 Test 4. Multiple sources in a cosmological 
density field 

The last test in RT06 involves a static cosmological density 
field at z = 9. The simulation comoving box size is 0.5 
Mpc and has a resolution of 128^. There are 16 point sources 
that are centred in the 16 most massive haloes. They emit 
fj = 250 ionising photons per baryon in a blackbody spec- 
trum with an effective temperature T = 10^ K, and they 
live for ts = 3 Myr. Thus the luminosity of each source is 



(37) 



where M is the halo mass, Qm = 0.27, and Qb — 0.043. The 
radiation boundaries are isolated so that the radiation leaves 
the box instead being shifted periodically. The simulation is 
evolved for 0.4 Myr. 

Fig. [ITl depicts the distribution of the neutral frac- 
tion and temperature of the entire domain and shows good 



5 RADIATION HYDRODYNAMICS TESTS 

We next show results from radia t ion hy drodynamics test 
problems presented in llliev et al.l (|2009l , hereafter RT09). 
They involve (1) the expansion of an H II region in a uniform 
medium, similar to Test 2, (2) an H II region in an isother- 
mal sphere, and (3) the photo-evaporation of a dense, cold 
clump, similar to Test 3. We turn off self- gravity and AMR 
in accordance with RT09. 



5.1 Test 5. Classical H ii region expansion 

Here we consider the expansion of an H ii region into a 
uniform neutral medium including the hydrodynamical re- 
sponse to the heated gas. The ionised region has a greater 
pressure than the ambient me dium, causing it to expand. 
This is a well-studied problem (|SDitzerlll978h with an ana- 
lytical solution, where the ionisation front moves as 
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Figure 20. Test 5. (H ll region in a uniform medium). Top: 
Growth of the computed ionisation front radius at an ionised 
fraction Xe = 0.5 (dashed) and at a temperature T = 10^ K (dot- 
ted) compared to the analytical estimate [solid; equation p8|) ]. 
Bottom: The ratio of the computed ionisation front radii to the 
analytical estimate. 



Figure 21. Test 5. (H ll region in a uniform medium). Clock- 
wise from the upper left: Radial profiles of density, temperature, 
ionised fraction, and pressure at times t = 10, 200, and 500 Myr. 



where Cs is the sound speed of the ionised gas and rs,o is the 
rs in Equation [35] The bubble eventually reaches pressure 
equilibrium with the ambient medium at a radius 

(39) 



4/7 



(38) 



/2T\2/3 
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where T and To are the ionised and ambient temperatures, 
respectively. These solutions only describe the evolution at 
late times, and not the fast transition from R-type to D-type 
at early times. 

The simulation setup is similar to Test 2 with the ex- 
ception of the domain size L = 15 kpc. Here pressure equi- 
librium occurs at r/ = 185 kpc, which is not captured by 
this test. However more interestingly, the transition from R- 
type to D-type is captured and occurs around i?s = 5.4 kpc. 
The test is run for 500 Myr. 

The growth of the ionisation front radius is shown in 
Fig. [20l using both T = 10^ K and Xe = 0.5 as ionisation 
front definitions, compared to the analytical solution [equa- 
tion (|39|) ]. We define this alternative measure because the 
ionisation front becomes broad as the D-type front creates 
a shock. Densities in this shock, as seen in Fig. I21[ are high 
enough for the gas to recombine but not radiatively cool. Be- 
fore 2trec ~ 250 Myr, the temperature cutoff overestimates 
Ts by over 10%; however at later times, it provides a good 
match to the t^^^ growth at late times. With the Xe = 0.5 
criterion for the ionisation front, the radius is always under- 
estimated by ^ 20%. This behavior was also seen in RT09. 

Fig. [21] shows the progression of the ionisation front at 
times t — 10, 200, and 500 Myr in radial profiles of den- 
sity, temperature, pressure, and ionised fraction. The ini- 
tial H II region is over-pressurised and creates an forward 
shock wave. The high-energy photons can penetrate through 
the shock and partially ionises and heats the exterior gas, 
clearly seen in the profiles. As noted in RT09, this heated 
exterior gas creates an photo-evaporative flow that flows in- 
ward. This interacts with the primary shock and creates the 
double-peaked features in the density profiles at 200 and 
500 Myr. Fig. 1221 shows slices through the origin of the same 
quantities, including neutral fraction. These depict the very 
good spherical symmetry of our method. The only appar- 
ent artifact is a very slight diagonal line, which is caused by 
the HEALPix pixelisation differences between the polar and 
equatorial regions. This artifact diminishes as the ray-to-cell 
sampling is increased. 

5.2 Test 6. H ii region expansion in an isothermal 
sphere 

A more physically motivated scenario is an isothermal 
sphere with a constant density Uc core, which is applica- 
ble to collapsing molecular clouds and cosmological haloes. 
The radial density profile is described by 



n(r) 



nc{r/ro) 



(r ^ ro) 
(r > ro) ' 



(40) 



where ro is the radius of the core. If the Stromgren radius is 
smaller than the core radius, then the resulting H ii region 
never escapes into the steep density slope. When the ioni- 
sation front propagates out of the core, it accelerates as it 
travels down the density gradient. There exists no analyti- 
cal solution for this pr oblem with fu l l gas d ynamics but was 
extensively studied bv [Franco et al.1 (|l990l l. After the gas is 
ionised and photo-heated, the density gradient provides the 
pressure imbalance to drive the gas outwards. 

This test is constructed to study the transition from 
R-type to D-type in the core and back to R-type in the den- 
sity gradient. Thus the simulation focuses on small-scale. 
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Figure 23. Test 6. (H ll region in a 1/r^ density profile). Top: 
Growth of the computed ionisation front radius as T = IC^ K 
and Xe = 0.5 as definitions for the front. Bottom: Velocity of 
the ionisation front, computed from outputs at 0.5 Myr intervals. 
The velocity is calculated from rip, whose coarse time resolution 
causes the noise seen in vip. It is smooth within the calculation 
itself. 
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Figure 24. Test 6. (H ll region in a 1/r^ density profile). Clock- 
wise from the upper left: Radial profiles of density, temperature, 
ionised fraction, and pressure at times t = 3, 10, and 25 Myr. 



not long term, behavior of the ionisation front. The simu- 
lation box has a side length L = 0.8 kpc with core density 
no = 3.2 cm~^, core radius ro = 91.5 pc (15 cells), and 
temperature T = 100 K throughout the box. The ionisa- 
tion fraction is initially zero, and the point source is located 
at the origin with a luminosity of 10^*^ ph s~^ cm~^. The 
simulation is run for 75 Myr. 

Because this problem does not have an analytical solu- 
tion, we compare our calculated ionisation front radius and 
velocity, shown in Fig. [23l to the RT09 results. Their evolu- 
tion are in agreement within 5% of RT09. As in Test 5, we 
use an extra definition of T = 10^ K for the ionisation front. 
We compute the ionisation front velocity from the radii at 
50 outputs, which causes the noise seen Fig. 1231 
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Figure 25. Test 6. (H ll region in a density profile). Clockwise from the upper left: Slices through the origin of ionised fraction, 

neutral fraction, temperature, and density at time t = 25 Myr. 



For the first Myr, the radiation source creates a weak R- 
type front where the medium is heated and ionised but does 
not expand because vw > Cs. When vw becomes subsonic, 
the medium can react to the passing ionisation front and 
creates a shock, leaving behind a heated rarefied medium. 
This behavior is clearly seen in the radial profiles of density, 
temperature, ionised fraction, and pressure in Fig. [24l The 
inner density decreases over two order of magnitude after 
25 Myr. To illustrate any deviations in spherical symmetry, 
we show in Fig. [25] slices of density, temperature, neutral 
fraction, and ionised fraction at the final time. The only ar- 
tifact apparent to us is the slight broadening of the shock 
near the x = and y — planes. This causes the ionisation 
front radius to be slightly smaller in those directions. In the 
diagonal direction, the neutral column density through the 
shock is slightly smaller, allowing the high-energy photons 
to photo-ionise and photo-heat the gas to Xe = 5 x 10~^ 
and T = 2000 K out to - 50 pc from the shock. The re- 
flecting boundaries are responsible for this artifact because 
this is not seen when the problem is centred in the domain, 
removing any boundary effects. 
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Figure 26. Test 7. (Photo-evaporation of a dense clump). Line 
cuts from the point source through the middle of the dense clump 
dit t = 1,10,50 Myr of (clockwise from the upper left) density, 
temperature, pressure, and neutral fraction. 



5.3 Test 7. Photo-evaporation of a dense clump 

The photo-evaporation of a dense clump in a uniform 
medium proceeds very differently when radiation hydrody- 
namics is considered instead of a static density field. The 



ionisation front first proceeds as a very fast R-type front, 
then it slows to a D-type front when it encounters the dense 
clump. As the clump is gradually photo- ionised and heated, 
it expands into the ambient medium. The test presented 
here is exactly like Test 3 but with gas dynamics. In this 
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Figure 27. Test 7. (Photo-evaporation of a dense clump). Clockwise from the upper left: Slices through the clump centre of neutral 
fraction, pressure, temperature, and density at time t = 10 Myr. 




Figure 29. Test 7. (Photo-evaporation of a dense clump). Prob- 
ability distribution function for temperature (left) and flow Mach 
number (right) at 10 Myr (solid) and 50 Myr (dashed). 



setup, the ionisation front overtakes the entire clump, which 
is then completely photo-evaporated. 

Fig.l26lshows cuts of density, temperature, neutral frac- 
tion, and pressure in a line connecting the source and the 
clump centre at t = 1, 10, and 50 Myr. At 1 Myr, the 
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Figure 30. Test 7. (Photo-evaporation of a dense clump). Slices 
through the clump centre of the flow Mach number at t = 10 Myr 
(left) and 50 Myr (right). 



ionisation front has propagated through the left-most 500 
pc of the clump. This heated gas is now over-pressurised, 
as seen in the pressure plot in Fig. [26l and then expands 
into the ambient medium. This expansion creates a photo- 
evaporative flow, seen in many star forming regions (e.g. 
M16; Hester et al .11 19961 1 as stars irradiate nearby cold, dense 
overdensities. These flows become evident in the density at 
10 Myr, seen both in the line cuts and slices fFig. [27)) . They 
have temperatures up to 50,000 K. At this time, the front has 
progressed about halfway through the clump, if one inspects 
the neutral fraction. However the high energy photons have 
heated all but the rear surface of the clump. At the end of the 
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Figure 28. Test 7. (Photo-evaporation of a dense clump). Same as Fig. [26]but at t = 50 Myr. 
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test (t = 50 Myr), only the core and its associated shadow 
is neutral, seen in Fig. [28l The core has been compressed by 
the surrounding warm medium, thus causing the higher den- 
sities seen at t = 50 Myr. The non-spherical artifacts on the 
inner boundary of the warm outermost shell are caused by 
the initial discretisation of the sphere, as discussed in ^4.31 
Next Fig.[29lshows the distributions of temperature and flow 
Mach number in the entire domain at 10 Myr and 50 Myr, 
showing similar behavior as the codes in RT09. Lastly in 
Fig. 1301 we show slices of the flow Mach number at 10 Myr 
and 50 Myr, showing the supersonic photo-evaporative flows 
that originate from the clump. 



6 RADIATION HYDRODYNAMICS 
APPLICATIONS 

We have completed presenting results from the RT06 and 
RT09 test suites. We expand on these test suites to include 
more complex situations, such as a Ray leigh- Taylor problem 
illuminated by a radiation source, champagne flows, an ir- 
radiated blast wave, collimated radiation, and an H ii with 
a variable source that further demonstrate its capabilities 
and accuracy. Last we use the new MHD implementation 
in ENZO v2.0 in the problem of a growing H ii region in a 
magnetic field. 



6.1 Application 1. 
clump 



Champagne flow from a dense 



Radiation-driven outflows from over densities , known as 
champagne flows, is a long studied problem fe.g .lYorkelll986l, 
§3.3). To study this, we use the same setup as lBisbas et al.l 
- a spherical tophat with an overdensity of 10 and 
radius of 1 pc in a simulation box of 8 pc. The ambient 
medium has p — 290 cm~^ and T — 100 K. The radiation 
source is offset from the overdensity centre by 0.4 pc. It 
has a luminosity of 10^^ ph s~^ and a T = 10^ K blackbody 
spectrum. The resulting Stromgren radius is 0.33 pc, just in- 
side of the overdense clump. These parameters are the same 
used in Bis bas et al. (2009). The entire domain initially has 
an ionised fraction of 10~^. We do not consider self-gravity. 
The simulation has a resolution of 128^ on base grid, and 
we refine the grid up to 4 times if a cell has an overdensity 
of 1.5 X 2^ where / is the AMR level. The simulation is run 
for 150 kyr. 

We show slices in the x-y and x-z planes of density in 
Fig. [SB at t = 10, 40, 100, 150 kyr. In the direction of the 
clump centre, the ionisation front shape transitions from 
spherical to parabolic after it escapes from the clump in 
the opposite direction. At t = 10 kyr, the surface of the H 
II region is just contained within the overdensity. In the x-z 
plane, there are density perturbations only above a latitude 
of 45 degrees. We believe that these are caused by the mis- 
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Figure 31. Application 1. (Champagne flow from a dense clump). Slices of density through the initial clump centre in the x-y plane (top) 
and x-z plane (bottom) at t = 10, 40, 100, 150 kyr. Notice the instabilities that grow from perturbations created while the H ll region is 
contained in the dense clump. 



match between HEALPix pixels and the Cartesian grid, even 
with our geometric correction. After the ionisation front es- 
capes from the clump in the negative x-direction, these per- 
turbations grow from Ray leigh- Taylor instabilities as the gas 
is accelerated when it exits the clump. As the shock propa- 
gates through the ambient medium, it is no longer acceler- 
ated and has a nearly constant velocity, as seen in Test 6. 
Thus these perturbations are not as vulnerable to Rayleigh- 
Taylor instabilities at this point. The ambient medium and 
shock ar e always optic ally thick, even in the directions of the 
bubbles. iBisbas et al. ^ found that the shock fragmented and 
formed globules; however we find the density shell is stable 
against such fragmentation. To investigate this scenario fur- 
ther, our next tests involve radiation driven Ray leigh- Taylor 
instabilities. 



6.2 Application 2. Irradiated Rayleigh- Taylor 
instability 

Here we combine the classic case of a Rayleigh- Taylor in- 
stability and an expanding H ii region. The Rayleigh- Taylor 
instability occurs when a dense fluid is being supported by 
a lighter fluid, initially in hydrostatic equilibrium, in the 
presence of a constant acceleration field. This classic test 
alone evaluates how subsonic perturbations evolve. We con- 
sider the case of a single-mode perturbation. The system 
evolves without any radiation until the perturbation grows 
considerably and then turn on the radiation source. These 
tests demonstrate that enzo+MORAY can follow a highly 
dynamic system and resolve fine density structures. 

We run two cases, an optically-thick and optically- thin 
case. In the fo rmer, w e take the parame ter choice s from past 
literature fe.g. lLiska Wendroff 2003: Stone et a l. 2008) by 
setting the top and bottom halves of the domain to a density 
yoi = 2 and po = 1, respectively. The velocity perturbation 



is set in the 2;-direction by 

Vz{x,y,z) = 0.01[1 + cos(27rx/La,)] X 
[1 + cos{27iy / Ly)] X 

[l + cos(27rz/L^)]/8. (41) 

We set the acceleration field Qz = 0.1 and the adia- 
batic index 7 = 1.4. We use a domain size of (^Lxt — 
(0.5,0.5, 1.5) with a resolution of (64, 64, 192). For hydro- 
static equilibrium, we set P = Po — gp{z)z with Po = 2.5. 
In order to consider a radiation source with a ionising pho- 
ton luminosity of 10^^ ph s~^, we scale the domain to a 
physical size of (0.5, 0.5, 1.5) pc; time is in units of Myr; 
density is in units of m^, resulting in an initial temperature 
of (To,Ti) = (363, 726) K. The radiation source is placed at 
the centre of lower 2;-boundary face and starts to shine at 
t = 10 Myr. 

The optically-thin case is set up similarly but with 
three changes: (1) a density contrast of 10, (2) a lumi- 
nosity of 10^^ ph s~^, and (3) the source is born at 6.5 
Myr. The time units are decreased to 200 kyr so that 
(To,Ti) = (1.8 X 10^, 1.8 X 10^) K. Note that in code units, 
pressure is unchanged. We adjust the physical unit scaling 
because we desire an optically thin bottom medium with 
T > 10^ K and Xe ^ 1. Furthermore, the ionisation front re- 
mains R-type before interacting with the instability. A pos- 
sible physical analogue could be a radiation source heating 
and rarefying the medium below. 

The X and ^/-boundaries are periodic, and the z- 
boundaries are reflecting. These will cause artificial features, 
in particular, because of the top reflecting boundary; never- 
theless, these tests provide a stress test on a radiation hy- 
drodynamics solver. We show the evolution of the density, 
temperature, and ionised fraction of the optically thick and 
optically thin cases in Figs. [32] and 1331 The initial state of 
the Rayleigh- Taylor instability is shown in the left panels. 

In the optically thick case, a D-type front is created. 
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Figure 32. Application 2. (Irradiated Ray leigh- Taylor instabil- 
ity; optically thick case). Slices at y = of density (top), temper- 
ature (middle), and electron fraction (bottom). The source turns 
on at t = 0. 



which is clearly illustrated by the spherical density enhance- 
ment at 0.02 Myr. The shock then passes through the insta- 
bility at ^0.25 Myr and reflects off the upper z-boundary. 
This and complex shock refl ections c reate a Richtmyer- 
Meshkov instability (see Brouillette 2002, for a review), driv- 
ing a chaotic jet-like structure downwards. The radiation 
source photo-evaporates the outer parts of this structure. 
The interaction between the dense cool "jet" and the hot 
medium further drives instabilities along the surface, which 
can be seen when comparing t = 0.59 Myr and t = 0.91 Myr 
slices. At the latter time, the jet cannot reach the bottom of 
the domain before being photo-evaporated. Eventually this 
structure is completely destroyed, leaving behind a turbu- 
lent medium between the hot and cold regions. 

The optically thin problem is less violent than the op- 
tically thick case because the R-type front does not interact 
with the initial instability as strongly. The radiation source 
provides further buoyancy in the already T = 10^ K gas. The 
gas first to be ionised and photo-evaporated are the outer 
regions of the instability. The enhanced heating also drives 
the upper regions of the instability, making the top interface 
turbulent. It then reflects off the upper z-boundary and cre- 
ates a warm T = 5 x 10^ K, partially ionised (xe ^ 10~^), 
turbulent medium, seen in the slices t ^ 0.67 Myr. The slices 
of electron fraction also show that the dense gas is optically 
thick. 




Figure 33. Same as Fig. [32]but for the optically thin case. 



6.3 Application 3. Photo-evaporation of a 
blastwave 

A supernova blastwave being irradiated by a nearby star is 
a likely occurrence in massive-star forming regions. In this 
test, we set up an idealised test that mimics this scenario. 
The ambient medium has a density po = 0.1 cm~^ and tem- 
perature To = 10 K. The domain size is 1 kpc. We use 2 levels 
of AMR with a base grid of 64^ that is refined if the density 
or total energy slope is greater than 0.4. The blastwave is 
initialised at the beginning of the Sedov- Taylor phase when 
the mass of the swept-up material equals the ejected mate- 
rial. It has a radius of 21.5 pc, a total energy of 10^° erg, and 
total mass of lOOM©, corresponding to ^ = 315 eV per par- 
ticle or E/kb — 3.66 x 10^ K. The radiation source is located 
at the centre of the left x-boundary and has a luminosity of 
10^° erg s"^ We use a T = 10^ K blackbody spectrum with 
2 energy groups (16.0 and 22.8 eV). The source turns on at 
2.5 Myr at which point the blastwave has a radius of 200 pc. 
The simulation is run for 7.5 Myr. 

Fig. [34] shows the ionisation front overtaking and dis- 
rupting the blastwave. We show the blastwave before the 
source is born at 2.5 Myr. The interior is rarefied [p ^ 
10~^ cm~^) and is heated to T ^ 5 x 10^ K by the reverse 
shock. At t = 3 Myr, the ionisation front is still R-type, 
and it ionises the rear side of the dense shell. Because the 
interior is ionised and diffuse, the ionisation front rapidly 
propagates through it until it reaches the opposite shell sur- 
face. Shortly afterward, the ionisation front transitions from 
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Figure 34. Application 3 (Photo-evaporation of a blastwave). Slices of density (top) and temperature (bottom) at t = 2.5,3,5,7.5 
Myr in the x — z plane. As the R-type ionisation front propagates through the blastwave centre, instabilities grow from the slightly 
inhomogeneous hot and rarefied medium. Note that the dense shell of the blastwave also creates dense inward fingers in the ionisation 
front shock. 
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Figure 35. Application 4 (Collimated radiation from a dense clump). Slices of density (top) and temperature (bottom) at t = 
0.1,3.25,10.75,23.25 Myr. The conical H ll region drives shocks transversely into the overdense sphere and creates polar champagne 
flows. The ambient medium is heated to T ~ 3 x 10^ K as the ionisation front passes the constant-pressure cloud surface. The ionisation 
front changes from D-type to R-type after it enters the ambient medium. 



R-type to D-type at a radius of 0.5 kpc, seen in the forma- 
tion of a shock in the 5 Myr density panel. This transition 
occurs by the construction of the problem not by the in- 
teraction with the blastwave. The surfaces of the blastwave 
that are perpendicular to the ionisation front have the high- 
est column density and thus are last to be fully ionised. The 
pressure forces from warm ambient medium and blastwave 
interior compress these surfaces, photo-evaporating them in 
the process, similar to Test 7. They survive until the final 
time t = 7.5 Myr. As the R-type ionisation front interacts 



with the blastwave interior, t he density perturbations c reate 
ionisation front instabilities (jWhalen Normanll2008h that 
are seen on the H ii region surface at the coordinate z — 0.5. 
Behind the ionisation front, the dense shell of the blastwave 
is photo-evaporated, and a smooth overdensity is left in the 
initial blastwave centre. 
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Figure 36. Application 4 (Time variations of the source luminos- 
ity). Top: Slice of the photo-ionisation rate fcph through the origin. 
The source has a duty cycle of 0.5 Myr, and the box has a light 
crossing time of 3.3 Myr. The shells of high kp]^ originate from 
radiation that was emitted when the source was at its peak lumi- 
nosity, illustrating the time-dependence of the radiative transfer 
equation. Bottom: Radial profile fcph with the inverse square law 
over plot ted. 



6.4 Application 4. Collimated radiation from a 
dense clump 

Some astrophysical systems produce collimated radiation ei- 
ther intrinsically by relativistic beaming or by an optically- 
thick torus absorbing radiation in the equatorial plane. The 
latter case would be applicable in a subgrid model of active 
galactic nuclei (AGN) or protostars, for example. Simulat- 
ing collimated radiation with ray tracing is trivially accom- 
plished by only initialising rays that are within some opening 
angle Oc- 

We use a domain that is 2 kpc wide and has an ambi- 
ent medium with po 10~^ cm"^, T 10^ K, Xe 0.99. 
We place a dense clump with p/po = 100, T = 100 K, 
Xe = 10~^, and r = 250 pc, at the centre of the box. Radi- 
ation is emitted in two polar cones with Oc = 7r/6 with 768 
(HEALPix level 3) initial rays, a total luminosity of 10^^ 
erg s~^, and a 17.6 eV mono-chromatic spectrum. This re- 
sults m tree — 

1.22 Myr and Rs — 315 pc, just outside of 
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Figure 37. Application 4 (Time variations of the source luminos- 
ity). Radial profiles of (clockwise from upper left) density, temper- 
ature, pressure, and ionised fraction at t = 50, 100, 150, 200, 250 
Myr. In this problem, the time variations in the source have little 
effect on the overall H ll region expansion. Inside the ionisation 
front at t = 50 Myr, there are small density perturbations that 
are created by the variable source that are later smoothed out 
over a sound crossing time. 

the sphere. The base grid has a resolution of 64^, and it is 
refined with the same overdensity criterion as Application 
1. We run this test for 25 Myr. 

We illustrate the expansion of the H ii region created 
by the beamed radiation in Fig. 1351 Before t = 3 Myr, the 
H II region is conical and contained within the dense clump, 
depicted in the t = 0.1 Myr snapshot of the system. At this 
time, the ionisation front is transitioning from R-type to D- 
type in the transverse direction of the cone. This can be seen 
in the minute overdensities on the H ii transverse surface. 
When it breaks out of the overdensity, a champagne flow de- 
velops, where the ionisation front transitions back to a weak 
R-type front. The cloud surface is a constant-pressure con- 
tact discontinuity (CD) with a density jump of 100. After 
the front heats the gas at the CD, there exists a pressure 
difference of ^ 100. In response, the high density gas accel- 
erates into the ambient medium and heats it to 3 x 10^ K. 
Additionally a rarefaction wave travels towards the clump 
centre. At later times, the transverse D-type front continues 
through the clump, eventually forming a disc-like structure 
at the final time. The polar champagne flows proceed to flow 
outwards and produces a dense shell with a diffuse (10~^^ 
cm~^) and warm (5000 K) medium in its wake. 

6.5 Application 5. Time variations of the source 
luminosity 

Our implementation retains the time derivative of the radia- 
tive transfer equation [equation (|2))] if we choose a constant 
ray tracing timestep, which saves the photon packages be- 
tween timesteps if c dtp < Lbox- This effect only becomes 
apparent when the variation time-scale of the point source is 
smaller than the light crossing time of the simulation. Fur- 
thermore, the timestep should resolve the variation time- 
scale by at least a few times. This property might be impor- 
tant in large box simulations with variable sources, e.g. AGN 



© 2010 RAS, MNRAS QQO.fWfl 



28 J. H. Wise and T. Abel 



20 



15 



:(pc) 
10 15 



20 



:(pc) 
10 15 



20 



3 10 



m 




\ 



5 • 10"^^ 

2 • 10-^^ 
1-10-22 1^ 

5 • 10-23 1: 

O 

2 • 10-23 ^ 

1 • 10-23 
5 • 10-24 



Figure 38. Application 5. (H ll region with MHD). Left to right: shces of density at t = 0.18,0.53,1.58 Myr in the x-y plane. The 
streamlines show the magnetic field. 
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Figure 39. Application 5. (H ll region with MHD). Slices of density (top) and the x-component of the magnetic field (bottom) in the 
y-z plane at t = 0.18,0.53, 1.58 Myr (left to right). 



radiative feedback. To test this, we can use an exponentially 
varying source with some duty cycle to - In a functional form, 
this can be described as 

L{t) = Lnaax X exp[A(t//to - 1)] (42) 

where t/ = 2x|t — toX round(t/to)|, and A = 4 con- 
trols the width of the radiation pulse. To illustrate the ef- 
fects of source variability, we remove any dependence on 
the medium by considering an optically-thin uniform den- 
sity p — 10~^ cm~^. We take Lbox — 1 Mpc, which has 
a light crossing time of 3.3 Myr. A source is placed at the 
origin with Lmax = 10^^ ph s~^ and to = 0.5 Myr. We use 
a radiative transfer timestep of 50 kyr to resolve the duty 



cycle by 10 timesteps. The simulation is run for 6 Myr, so 
the radiation propagates throughout the box. 

The variability of the source is clearly illustrated in the 
photo-ionisation rates, shown in Fig. 1361 The shells of rela- 
tive maximum kph corresponds to radiation that was emit- 
ted when the source was at its peak luminosity. They are 
separated by cto Mpc and are geometrically diluted with 
increasing radius. Averaged over shells of the same width, 
photo-ionisation rates decrease as 

Next we test the hydro dynamical response to a vary- 
ing source by repeating Test 5. We set the peak luminosity 
Lma^ = 2 X 10^^ erg s~^ that is a factor of 4 more luminous 
than Test 5, so the average luminosity is ^ 5 x lO'^^ erg s~^. 
The spectrum is mono-chromatic with an energy of 29.6 eV. 
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We set the variation time-scale t/ = cLbox/3 = 16.3 kyr and 
use a constant radiative transfer timestep tp = t//4 = 4.07 
kyr. The simulation is run for 200 Myr. We show the ra- 
dial profiles of density, temperature, ionised fraction, and 
pressure in Fig. 1371 The variable source has little effect on 
the overall growth of the H ii region. It has the approxi- 
mately the same radius as Test 5 at t = 200 Myr when run 
with a mono-chromatic spectrum (see ^7.3p . At early times, 
the variable source creates density perturbations with an 
average size of 500 pc inside the ionised region, seen in the 
t = 50 Myr profiles. They do not create any instabilities and 
are smoothed out over its sound crossing time of ^ 50 Myr. 



6.6 Application 6. H ii region with MHD 

Another prevalent physical component in astrophysics is a 
magnetic field. We util ise the new ma^net ohydrodynam- 
ics (MHD) framework ("Wan g k Abel 12009') in ENZO v2.0 
that uses an unsplit conservative hydrodynamic s solver and 
the hyperbolic V • B = cleaning method of iDedner et al.l 
([2002). We show a test problem with an expanding H 
II region in an initially uniform density field and con- 
stant magnetic f ield. We use the same problem setup as 
iKrumholz eraJ] (|2007l l - p = 100 cm"^ T = 11 K, Lbox 
20 pc with a resolution of 256^. This ambient medium is 
threaded by a magnetic field B = 14. 2x fiG. The Alfven 
speed is 2.6 km s~^. The radiation source is located in 
the centre of the box with a luminosity L = 4 x 10^^ 
ph s~^ with a 17.6 eV mono-chromatic spectrum, result- 
ing in a Stromgren radius Rs = 0.5 pc. The simulation is 
run for 1.58 Myr . The hydrodynami cs solver uses an HLL 
Riemann solver (|Harten et al.l 1 19831) and piec ewise linear 
method (PLM) reconstruction (va n Lee'rlll977l ) for the left 
and right states in this problem. 

As the H II region grows in the magnetised medium, 
shown in Figs. [38] and 1391 it transforms from spherical to 
oblate as it is magnetically confined in directions perpendic- 
ular to the magnetic field. This occurs at t > 0.5 Myr be- 
cause the magnetic pressure exceeds the t hermal pressure, 
and the gas can only flow along field lines. iKrumholz et al I 
observed some carbuncle artifacts along the ionisation front; 
whereas we see smooth density gradients, which is most 
likely caused by both the geometric correction to the ray 
tracing ( ^2.3p and the diffusivity of the HLL Riemann 
solver when conipare d to Roe's Riemann solver used in 
IKrumholz et al.l (|2007 ). who also use PLM as a reconstruc- 
tion method. The evolution of the magnetic field lines evolve 
in a similar manner as their results. 



7 RESOLUTION TESTS 

Resolution tests are important in validating the accuracy 
of the code in most circumstances, especially in produc- 
tion simulations where the initial environments surround- 
ing radiation sources are unpredictable. In this section, we 
show how our adaptive ray-tracing implementation behaves 
when varying spatial, angular, frequency, and temporal res- 
olutions. 









-' 1 ; 
/ 


16^ 

- - - 32^ 

64^ 

128^ 



100 200 300 400 500 

Time (Myr) 



Figure 40. Growth of the ionisation front radius, compared to 
the analytical radius, in Test 1 with varying spatial resolutions. At 
resolutions of 16^ and 32^, the ionisation front is underestimated 
for the first ~ 25 Myr but converges within 0.5% of the higher 
resolution runs. 

7.1 Spatial resolution 

Here we use Test 1 f ^.ip as a testbed to investigate how the 
evolution of the Stromgren radius changes with resolution. 
We keep all aspects of the test the same, but use resolutions 
of 16^ 32^ 64^ and 128^ In Fig. |40l we show the ratio 
^ip/Tanyi, similar to Fig. O using these different resolutions. 
The radii in the 64^ and 128^ runs evolve almost identi- 
cally. Compared to these resolutions, the lower 16^ and 32^ 
resolution runs only lag behind by 1% until 300 Myr, and 
afterwards it is larger by 0.5% than the higher resolution 
cases. This shows that our method gives accurate results, 
even in marginally resolved cases, which is expected with a 
photon conserving method. Furthermore this demonstrates 
that the geometric correction does not significantly affect 
photon conservation. 

7.2 Angular resolution 

The Cartesian grid must been sampled with sufficient rays in 
order to calculate a smooth radiation field. To determine the 
dependence on angular resolution, we consider the propaga- 
tion of radiation through an optically thin, uniform medium. 
The radiation field should follow a 1/r^ profile. As the grid 
is less sampled by rays, the deviation from 1/r^ should in- 
crease. This test is similar to Test 1, but the medium has 
p = 10"^ cm"^, T = 10"^ K, and 1 - Xe = 10"^. The simula- 
tion is only run for one timestep because the radiation field 
should be static in this optically-thin test. 

We consider minimum ray-to-cell ratios $c = 
(1.1,2.1,3.1,5.1,10.1,25.1). Slices of the photo-ionisation 
rates through the origin are shown in Fig. Hi] for these val- 
ues of $c. In this figure, we limit the colormap range to a 
factor of 3 to show the nature of the artifacts in more con- 
trast. Unsealed, the rates in the figures would span 4 orders 
of magnitude. When <I>c ^ 3.1, the cell-to-cell variations are 
apparent because there are not enough rays to sufficiently 
sample the radiation field, even with the geometric correc- 
tion factor /c, whose improvements are shown later in ^8.11 
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Figure 41. Variations in the photo-ionisation rates for different 
ray-to-cell samplings ^c- The color map only spans a factor of 3 
to enhance the contrast. In comparison, the photo-ionisation rate 
actually spans 4 orders of magnitude in this test. 
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Figure 42. Standard deviations of the difference between the 
computed photo-ionisation rates and an inverse square law as a 
function of ray-to-cell samplings <l>c for different spatial resolu- 
tions. There is no dependence on the spatial resolution, and the 
accuracy increases as a oc 



At $c = 5.1, these artifacts disappear, leaving behind a shell 
artifact where the radiation fields do not smoothly decrease 
as 1/r^. At higher values of $c, this shell artifact vanishes 
as well. 

One measure of accuracy is the deviation from an 1/r^ 
field because this problem is optically-thin. To depict the 
increase in accuracy with ray sampling, we take the differ- 
ence between the calculated photo-ionisation rate and a 1/r^ 
field, and then plot the standard deviation of this difference 
field versus angular resolution in Fig. [42l We plot this rela- 
tion for resolutions of 32"^, 64^, and 128^ and find no depen- 
dence on spatial resolution, which is expected because we 
control the angular resolution in terms of cell widths, not 
in absolute solid angles. We find that the deviation from an 
inverse square law decreases as a oc 
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Figure 43. Radial profiles of (clockwise from the upper left) 
density, temperature, radial velocity, and ionised fraction for Test 
5 with rijy = 1, 2, 4, 8, and 16 frequency bins sampling the T = 10^ 
K blackbody spectrum. The data are shown at t = 200 Myr (top) 
and t = 500 Myr (bottom). The double-peaked structure in the 
shock only appears with a multi-frequency spectrum. The solution 
converges at ^ 4. 



7.3 Frequency resolution 

The ionisation front radius is within 5-10% of analytical 
solutions in Tests 1, 2, and 5 with only one energy group; 
however a multi-frequency spectrum can create differences 
in the reactive flows. We use Test 5 f ^5.1[ an expanding H 
II region with hydrodynamics) to probe any differences in 
the solution when varying the resolution of the spectrum. 
In RT09, ZEUS-MP was used to demonstrate the effect of 
a multi-frequency spectrum on the dynamics of the ioni- 
sation front in this test. Instead of a single shock seen in 
the mono-chromatic spectrum, the shock obtains a double- 
peaked structure in density and radial velocity. We rerun 
Test 5 with a T = 10^ K blackbody spectrum sampled by 
riu = 1, 2, 4, 8, and 16 frequency bins. We use the following 
energies: 

• riiy = 1: Mean energy of 29.6 eV 

• riiy — 2: Mean energies in bins 13.6-30 and >30 eV-21.1, 
43.0 eV 
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Figure 44. Growth of the ionisation front radius, compared to 
the analytical radius, in Test 1 with varying radiative transfer 
timesteps. The dnn/dt and dl/dt based timesteps provide the 
best accuracy, combined with computational efficiency because 
they take short timesteps when the H ll is expanding rapidly but 
take long timesteps when the photon gradients are small when 
riF is large. At the final time, all but the t = 5 Myr constant 
timestep produce identical ionisation front radii. 



• = 4: Mean energies in bins 13.6-20, 20-30, 30-40, 
and >40 eV-16.7, 24.6, 34.5, 52.1 eV 

• riu = 8, 16: Logarithmically spaced between 13.6 and 
50 eV for the first riu — 1 bins, and the last bin is the mean 
energy above 50 eV. 

Fig- US] shows the radial profiles of density, temperature, 
ionised fraction, and radial velocity at t = 200 Myr and 
t = 500 Myr. All of the runs with > 1 show the double 
peaked features in density and radial velocities. The mono- 
chromatic spectrum misses this feature completely because 
all of the radiation is absorbed at a characteristic column 
density. In the multi-frequency spectra, the higher energy 
photons are absorbed at larger column densities and photo- 
heated this gas. This heated gas creates a photo-evaporative 
flow that collides with the innermost shock, forming the 
double peaked density profile. The n^y ^ 4 runs are in- 
distinguishable, and the riu = 2 spectrum only leads to 
a marginally higher density in the outer shock and lower 
ionised fractions and temperatures in the ambient medium. 
In effect, a mono-chromatic spectrum can be sufficient if the 
problem focuses on large-scale quantities, e.g. ionised filling 
fractions in reionisation calculations. Conversely these ef- 
fects may be important when studying the details of small- 
scale processes, e.g. photo-evaporation. 

7.4 Temporal resolution 

The previous three dependencies did not affect the propaga- 
tion of the ionisation fro nt greatly. However in our and oth- 
ers' past experience fe.g.[ShaDiro et al.ll200"3 : iMellema et al.l 
120061 : IPetkova SDringelll2009l \ the timestep, especially too 
small a one, can drastically underestimate the ionisation 
front velocity. Here we use Test 1 but with 64^ resolu- 
tion to compare different time-stepping methods - restricted 
changes in H ii (dnu/dt based; ^3.4.ip , constant timesteps of 



Figure 45. Variable time-stepping for the methods that limit 
change in neutral fraction (solid) and specific intensity (dotted). 
The horizontal lines show the constant timesteps that were used in 
the tests. The crossing time of a mean free path by the ionisation 
front is plotted for reference. 

0.1, 0.5, 1, and 5 Myr f ^3.4.3p , and based on incident radia- 
tion {dl/dt based; ^3.4.4p . The growth of the ionisation front 
radius is shown in Fig. 1441 Both the H ii restricted and inci- 
dent radiation variable time-stepping methods agree within 
a few percent throughout the entire simulation, as does the 
run with constant dtp = 0.1 Myr timesteps. With the larger 
constant timesteps, the numerical solution lags behind the 
analytical one, but they converge to an accurate H ii radii 
at late times. Even dtp = 5 Myr timestep, which under- 
estimated it by 35% at 50 Myr, is within a percent of the 
analytical solution. 

The larger constant timesteps deviate from these more 
accurate solutions at early times because the photon en- 
ergy gradient is large, and thus so is the ionisation front 
velocity. To understand this, the ionisation front can be 
considered static in a given timestep. Here the ionising ra- 
diation can only penetrate into the neutral gas by roughly 
a photon mean path Amfp- Only in the next timestep, the 
ionisation front can advance. If the timestep is larger than 
Amfp/'^^iF,anyi, thcu the uumcrical solution may fall behind. 

The variable time-stepping of the dnu/dt and dl/dt 
methods adjust accordingly to the physical situation, as seen 
in the plot of timestep versus time (Fig. |45)). They pro- 
vide high accuracy when the source first starts to shine. At 
later times, the ionisation front slows as it approaches the 
Stromgren radius, and large timesteps are no longer neces- 
sary. The dl/dt method has a similar timestep as the dnn/dt 
method. It is larger by a factor of ^ 2 because of our choice 
of the safety factor CRT,cfl = 0.5. This causes its calculated 
radius to be smaller by 1% at t < tree, which is still in good 
agreement with the analytical value. 



8 METHODOLOGY TESTS 

Here we show tests that evaluate new features in 
ENZO+MORAY, such as the improvements from the geometric 
correction factor, optically- thin approximations, treatment 
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Figure 46. Slices of the photo-ionisation rate in the x-y plane 
(top row) and x-z plane (bottom row) with (left column) and 
without (right column) the geometric correction. The slices are 
through the origin. In the x-y plane, it reduces the shell artifacts. 
In the x-z plane, it reduces the severity of a non-spherical arti- 
fact delineated at a 45 degree angle, where the HEALPix scheme 
switches from polar to equatorial type pixels. 

of X-ray radiation, and radiation pressure. Lastly we test for 
any non-spherical artifacts in the case of two sources. 




Figure 47. Optically-thin approximation to the radiation field 
with one ray per cell in optically thin regions. The angular arti- 
facts result from the transition to optically thick (white line) at 
an optical depth r = 0.1. 

upper region, they are polar HEALPix pixels. This artifact 
is not seen in the x-y plane because all rays are of a equa- 
torial type. The geometric correction smooths this artifact 
but does not completely remove it. 



8.1 Improvements from the covering factor 
correction 

As discussed in ^2.31 non-spherical artifacts are created by a 
mismatch between the HEALPix pixelisation and the Carte- 
sian grid. This is especially apparent in optically-thin re- 
gions, where the area of the pixel is greater than the (1 — e~^) 
absorption factor. In this section, we repeat the angular res- 
olution tests in ^7.21 with $c = 5.1. Slices of the photo- 
ionisation rates through the origin are shown in Fig. I46[ de- 
picting the improvements in spherical symmetry and a closer 
agreement to a smooth 1/r^ profile. Previous attempts to 
reduce these artifacts either introd uced a random rotation 
of the HEALPix pixelisation (e.g. lAbel Wandeltl liooi : 
iTrac CenI l2007l : iKnimholz et all l2007l ) or by increasing 
the ray-to-cell sampling. 

In the x-y plane without the correction, there exists 
shell artifacts where the photo-ionisation rates abruptly 
drops when the rays are split. This occurs because the pho- 
ton flux in the rays are constant, so kph is purely dependent 
on the ray segment length through each cell. Geometric dilu- 
tion mainly occurs when the number of rays passing through 
a cell decreases. With the correction, geometric dilution also 
occurs when the ray's solid angle only partially covers the 
cell. This by itself alleviates these shell artifacts. In the x-z 
plane without the correction, there is a non-spherical artifact 
delineated at a 45 degree angle. In the lower region, the rays 
are associated with equatorial HEALPix pixels, and in the 



8.2 Optically-thin approximation 

In practice, we have found it difficult to transition from the 
optically-thin approximation to the optically-thick regime 
without producing artifacts in the photo-ionisation rate kph- 
We use the optically thin problem used in the angular res- 
olution test f ^7.2|) with $c = 5.1 to show these artifacts in 
kph in Fig. [471 The radiation field strictly follows a 1/r^ pro- 
file until it reaches Tthin = 0.1, which is denoted by the white 
quarter circle in the figure. Within this radius, only $c = 1 
is required. Then at this radius, the rays are then split until 
a sampling of $c is satisfied. Angular spike artifacts be- 
yond this radius arise because of the interface between the 
optically-thin approximation and full ray tracing treatment. 
They originate in cells that intersect the Tthin surface, which 
are split into the optically thin and thick definitions. Unfor- 
tunately we have not determined a good technique to avoid 
such artifacts. They occur because of the following reason. 
When the first ray with r < 0.1 exits such a cell, it applies 
the optically-thin approximation and marks the cell so no 
other ray from the same source contributes to its kph and F 
field. However other rays may exit the cell with r > 0.1 be- 
cause the maximum distance between the far cell faces and 
the source is not always r < 0.1. Then these rays will split in 
this cell and add to kph and become attenuated, reducing its 
photon fiuxes. When the rays continue to the next cell after 
this transition, the photon fiuxes are not necessarily equal to 
each other, creating the angular artifacts seen in Fig.[47l We 
are continuing to formulate a scheme that avoids these arti- 
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Figure 48. Radial profiles of temperature and ionised frac- 
tion showing the effects of secondary ionisations from a mono- 
chromatic 1 keV spectrum. The discontinuity at r ~ 0.4 is caused 
by artifacts in the ray tracing, which is described in Test 1. The 
high energy photons can ionise multiple hydrogen atoms, increas- 
ing the ionised fraction. In part, less radiation goes into thermal 
energy, lowering the temperature. 

facts because this approximation will be very advantageous 
in simulations with large ionised filling factors. 

8.3 X-Ray secondary ionisations and reduced 
photo-heating 

Here we test our implementation of secondary ionisations 
from high-energy photons above 100 eV, described in ^2.5.21 
and used in Alva rez et al. (2009) in the context of accreting 
black holes. We use the same setup as Test 5 but with an in- 
creased luminosity L = 10^° erg s~^ and a mono-chromatic 
spectrum of 1 keV. Fig. [48] compares the density, tempera- 
ture, ionised fraction, and neutral fraction of the expanding 
H II region considering secondary ionisations and reduced 
photo-heating and considering only one ionisation per pho- 
ton and the remaining energy being thermalised. 

Fig. [48] shows the main effects of secondary ionisations 
from the 1 keV spectrum on the ionisation and thermal 
state of the system. Without secondary ionisations, each ab- 
sorption results in one ionisation with the remaining energy 
transferred into thermal energy. But with secondary ioni- 
sations, recall that most of the radiation energy goes into 
hydrogen and helium ionisations in neutral gas; whereas in 
ionised gas, most of the energy is thermalised. In this test, 
only the inner 300 pc is completely ionised because of the 
small cross-section of hydrogen at = 1 keV. Beyond this 
core, the medium is only partially ionised. This process ex- 
pands the hot T = 10^ K core by a factor of 2. In the outer 
neutral regions, the ionisation fraction is larger by a factor 
of ~ 10, which in turn results in less photo- heating, lowering 
the temperature by a factor 2-3. 

8.4 Radiation pressure 

Radiation pressure affects gas dynamics in an H ii region 
when its force is comparable to the acceleration created by 




^/^box ^/^box 

Figure 49. (a) No radiation pressure. Radial profiles of (clock- 
wise from top left) density, temperature, radial velocity, and neu- 
tral fraction. Time units in the legend are in kyr. (b) Radial pro- 
files with radiation pressure. The momentum transferred to the 
gas drives out the gas at higher velocities than without radiation 
pressure. Afterwards the central region is under-pressurised, and 
the gas infalls toward the centre, as seen at t = 60 kyr. Then the 
radiation pressure continues to force the gas outwards, increasing 
gas velocities up to 50 km/s. 



gas pressure of the heated region. The imparted accelera- 
tion on a hydrogen atom arp = ^ph/c. This is especially 
important when the ionisation front is in its initial R-type 
phase, where the gas has not reacted to the thermal pres- 
sure yet. Thus we construct a test that focuses on a small 
scales, compared to the Stromgren radius. The domain has 
a size of 8 pc with a uniform density p = 2900 cm~^ and 
initial temperature T — 10^ K. The source is located at the 
origin with a luminosity L — 10^° ph s~^ and a T = 10^ K 
blackbody spectrum. We use one energy group ^ph = 29.6 
eV. The grid is adapt ively refined on overdensity with the 
same criterion as Test 8. The simulation is run for 140 kyr. 

We compare nearly identical simulations but one with 
radiation pressure and one without radiation pressure to 
quantify its effects. Radial profiles of density, temperature, 
neutral fraction, and radial velocity are shown in Fig. 1491 for 
both simulations at several times. Without radiation pres- 
sure, the evolution of the H ii is matches the analytical 
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Figure 50. Slice of neutral fraction at t = 500 Myr through 
the sources in the consolidated H ll test. There are no artifacts 
associated with rays being emitted from two sources. Both of the 
ionised regions are spherically symmetric before they overlap. 



expectations described in ^5.11 At t = 140 kyr with radi- 
ation pressure, the ionisation front radius is increased by 
^ 5% = 0.16 pc. However radiation pressure impacts the 
system the greatest inside the ionisation front. At t = 40 
kyr, the central density is smaller by a factor of 20 with 
radiation pressure, but the temperatures are almost equal. 
A rarefaction wave thus propagates toward the centre, de- 
picted by the negative radial velocities at t = 60 kyr. This 
raises the central density to 10 g cm~^ at t = 80 kyr. 
Afterwards, the radiation continues to force gas outwards. 
From t — 100 kyr to t = 140 kyr, the maximum radial veloc- 
ity of the ionised gas increases from 10 km s~^ to 50 km s~^. 
This leaves behind an even more diffuse medium, lowering 
the gas density by a factor of 10 at t = 140 kyr. Thus the re- 
combination rates are lower, resulting in increased ionisation 
fractions and temperatures in the H ii region. 



8.5 Consolidated H ii region with two sources 

Here we test for any inaccuracies in the case of 
multiple sources. We u se the same test problem as 
IPetkova Springel (|2009l , §5.1.2), which has two sources 
with luminosities of 5 x 10^^ ph s~^ and are separated by 8 
kpc. The ambient medium is static with a uniform density 
of 10"^ cm"^and T = 10^ K. This setup is similar to Test 
1. The domain has a resolution of 128 x 64 x 64 It is 20 kpc 
in width and is 10 kpc in height and depth. The problem is 
run for 500 Myr. 

The H II regions grow to r = 4 kpc where they overlap. 
Then the two sources are enveloped in a common, elongated 
H II region. To illustrate this, we show the neutral fraction 
in Fig. 1501 Our method keeps spherical symmetry close to 
the individual sources, and there are no perceptible artifacts 
from having multiple sources. 



9 PARALLEL PERFORMANCE 

Last we demonstrate the parallel performance of 
ENZO+MORAY in weak and strong scaling tests. For large 
simulations to consider radiative transfer, it is imperative 
that the code scales to large number of processors. 



9.1 Weak Scaling 

Weak scaling tests demonstrates how the code scales with 
the number of processors with a constant amount of work 
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Figure 51. Weak scaling test with one 64^ block per process. 
Each block has one source and is set up similar to the radiation 
hydrodynamics Test 5. Above 8 processes, all parts of the code 
exhibit good weak scaling except for the inter-processor ray com- 
munication. The radiation module timing include the ray tracing, 
communication, chemistry and energy solver, and all other over- 
heads associated with the radiation transport. Cache locality of 
the data causes the decrease in performance from 1 to 8 processes. 



per processor. Here we construct a test problem with a 64^ 
block per core. The grid is not adapt ively refined. The phys- 
ical setup of the problem is nearly the same as Test 5 with 
a uniform density p = 10~^ cm~^ and initial temperature 
T = 100 K. Each block has the same size of 15 kpc as Test 
5. At the centre of each grid, there exists a radiation source 
with a luminosity L = 5 x 10^^ ph s~^ and a 17 eV mono- 
chromatic spectrum. The problem is run for 250 Myr. We 
run this test with Np = 2^ cores with n = [0, 1 ... 10, 11]. 
The domain has {Nx, Ny, Nz) blocks that is determined with 
the MPI routine MPI_Dims_create. For example with n — 7, 
the problem is decomposed into {Nx, Ny, Nz) = (4,4,8) 
blocks, producing a 256 x 256 x 512 grid. We have run these 
on the original (Harpertown CPUs) nodes of the NASA NAS 
machine, Pleiades, with 8 cores per node. Fig. [511 shows the 
performance timings of various parts of the code. From one 
to two cores, the total time only increases by ^1% due to the 
overhead associated with the inter-processor message pass- 
ing. From two to eight cores, the performance decreases in 
the hydrodynamics solver, chemistry and energy solver, and 
obtaining the hydrodynamic boundary conditions because 
of cache locality problems of the data being passed to the 
CPU. This occurs because of the CPU architecture, specif- 
ically the Ll cache and core connectivity. We see less of a 
penalty in newer processors, e.g. Intel Nehalam and West- 
mere CPUs. Above eight processes, these routines exhibit 
near perfect weak scaling to 2048 cores. Unfortunately there 
exists a Np^ dependence in the ray communication routines. 
It becomes the dominant process above 512 processes. We 
are actively pursing a solution to this scaling problem. The 
other parts of the code exhibit excellent weak scaling. Over- 
all it scales already well to 512 processes, and there is room 
to enhance the weak scalability of enzo+moray in the near 
future. 
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Figure 52. Strong scaling test with a 256^ AMR calculation of 
cosmological reionisation with ~ 392^ zones. The error bars rep- 
resent the minimum and maximum time spent on a single core 
and measures load balancing. Each point has been slightly off- 
set to make the error bars distinguishable. The hydrodynamics 
and non-equilibrium chemistry solvers scale well. The communi- 
cation of rays does not scale well in this problem and limits the 
performance, where as the ray tracing scales relatively well. 



9.2 Strong Scaling 

Strong scaling tests shows how the problem scales with the 
number of processors for the same problem. The overhead 
associated with the structured AMR framework in enzo can 
limit the strong scalability. One key property of strong scal- 
ing is that each processor must have sufficient work to com- 
pute, compared to the communication involved. In our ex- 
perience, non-AMR calculations exhibit much better strong 
scaling than AMR ones because of reduced inter-processor 
communication. We use an AMR simulation to demonstrate 
the scalability of enzo+moray in a demanding, research ap- 
plication. Here we use a small-box reionisation calculation 
with Lbox = 3 Mpc/h, a resolution of 256^, and six levels 
of refinement. We measure the time spent on the hydrody- 
namics, non-equilibrium chemistry, ray tracing, and radia- 
tion transport communication in a timestep lasting 1 Myr, 
at z = 10. There are nine radiative transfer timesteps in this 
period. The box has 675 point sources, 15,943 AMR grids, 
and 6.01 x 10^ ~ 392^ computational cells in this calcula- 
tion at this redshift. On average, 4.6 x 10^ ray segments are 
traced each timestep. The ionised volume fraction is 0.10. 
The calculations are performed on the Nehalam nodes on 
Pleiades on 2^ cores, where n = [2, 9]. 

Fig. [52] shows the strong scaling results of this calcu- 
lation. The hydrodynamics and non-equilibrium chemistry 
routines scale very well in this range because they depend 
on local phenomena. Note that the dominant process is the 
radiation transport instead of the hydrodynamics when com- 
pared to the weak scaling tests. The error bars in the figure 
represent the deviations across all cores, and it shows that 
load balancing becomes an issue at 128 and 256 cores. This 
is an even larger problem for the ray tracing because the 
grids that host multiple radiation sources will have signifi- 
cantly more work than a distant grid. For example if a few 
central grids host several point sources, the exterior grids 



must wait until those grids are finished ray tracing in order 
to receive the rays. This idle time constitutes the majority 
of the time spent in the ray communication routines despite 
our efforts (see ^3.5p to minimize idle time. For this reason, 
the communication of rays does not scale in this problem. 
One solution is to split the AMR grids into smaller blocks 
based on the ray tracing work. This could impact the per- 
formance of the rest of enzo by increasing the number of 
boundaries and thus communication. Nevertheless in sim- 
ulations that are dominated by radiation transport, it will 
be advantageous to construct such a scheme to increase the 
feasibility of running larger problems. 



10 SUMMARY 

In this paper, we have presented our implementation, 
ENZO+MORAY, of adaptive ray tracing ( Abel Wan deltl 
2002) and its coupling to the hydrodynamics in the cosmol- 
ogy AMR code enzo, making it a fully functional radiation 
hydrodynamics code. As this method is photon conserving, 
accurate solutions are possible with coarse spatial resolu- 
tion. A new geometric correction factor to ray tracing on 
a Cartesian grid was described, and it is general to any 
implementation. We have exhaustively tested the code to 
problems with known analytical solutions and the problems 
presented in the RT06 and RT09 radiative transfer compari- 
son papers. Additionally we have tested our code with more 
dynamical problems - champagne flows, Ray leigh- Taylor in- 
stabilities, photo-evaporation of a blastwave, beamed radia- 
tion, a time- varying source, and an H ii region with MHD - 
to demonstrate the flexibility and fidelity of enzo+moray. 
Because production simulations may not have the resolu- 
tion afforded in these test problems, we have tested the 
dependence on spatial, angular, frequency, and temporal 
resolution. It provides accurate solutions even at low res- 
olution, except for the large constant timesteps. However, 
we have described two methods to determine the radiative 
transfer timestep that are based on the variations in specific 
intensity or changes in neutral fraction inside the ionisation 
front. Both methods give very accurate results and provide 
the largest timestep to obtain an accurate solution, ulti- 
mately leading to higher computational efficiency. On the 
same topic, we have described a method to calculate the 
radiation field in the optically-thin limit with ray tracing. 

Being a ray tracing code, it scales with the number of 
radiation sources; nevertheless, it scales well to 0(10*^) pro- 
cessors for problems with ~ 10^ computational cells and 
^ lO'^ sources, such as reionisation calculations. We have 
also shown that the code shows good strong scaling in AMR 
calculations, given a large enough problem. The combination 
of AMR and adaptive ray tracing allows for high-resolution 
and high-dynamical range problems, e.g. present-day star 
formation, molecular cloud resolving cosmological galaxy 
formation, and H ii regions of Population III stars. Fur- 
thermore, we have included Lyman- Werner absorption, sec- 
ondary ionisations from X-ray radiation, Compton heating 
from photon scattering, and radiation pressure into the code, 
which extends the reach of enzo+moray to study AGN 
feedback, stellar winds, and local star formation. Coupling 
the radiative transfer with MHD further broadens the appli- 
cability of our code. The full implementation is included in 
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the latest public version of ENZclfl, providing the community 
with a full- featured radiation hydrodynamics AMR code. 
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